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CHAPTER  I 
INTRODUCTION 


Background 

Seismic  refraction  experiments  have  been  used  extensively  in  the  past 
thirty  five  years  in  investigations  of  the  structure  of  the  oceanic  crust.  The 
longer  range  of  the  refraction  or  wide  angle  reflection  technique,  on  the 
order  of  tens  of  kilometers,  permits  a  deeper  and  wider  area  of  examination, 
although  with  less  resolution,  than  the  spatially  limited  seismic  reflection 
experiment.  Observations  of  arrivals  from  the  Mohorovicic  discontinuity,  at  an 
average  depth  of  seven  kilometers  below  the  sea  floor,  are  routinely  made. 

The  major  focus  in  interpreting  refraction  data  has  been  the  analysis  of 
travel  time/range  data  and  the  "inversion"  of  this  data  for  the  purpose  of 
determining  a  velocity  versus  depth  profile  of  the  crust.  The  most  frequent 
application  of  this  procedure  is  the  geophysicist's  use  of  velocities  for 
postulating  geologic  structures  and  rock  types  below  the  sediment  (Christensen 
&  Salisbury,  1975).  Another  area  using  refraction  data,  less  widely  seen, 
falls  into  the  ocean  acoustician's  domain.  In  studying  the  behaviour  of  sound 
in  the  ocean,  the  sea  floor  is  often  modelled  as  a  boundary  with  a  half  space 
below,  and  with  some  form  of  reflection  characteristic  and/or  loss  mechanism. 
If  acoustic  energy,  upon  encountering  the  bottom,  was  either  reflected  or 
transmitted  directly,  this  would  be  appropriate,  and  the  determination  of 
reflection  and  transmission  coefficients  for  the  sea-sediment  interface  would 
probably  be  sufficient.  However,  sound  energy  does  penetrate  beneath  the  sea 
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floor  and  is  Doth  reflected  and  refracted  back  to  the  water.  In  an  active 
acoustical  experiment,  especially  at  longer  ranges,  a  significant  amount  of 
the  received  energy  may  come  from  waves  that  have  interacted  with  the  earth's 
crust  and  have  been  reinjected  into  the  water.  Since  these  arrivals  can  be 
detected  in  the  ocean,  their  study  is  of  concern  for  the  acoustician. 

The  role  of  bottom  interaction,  especially  at  low  frequencies,  is  now  an 
area  of  intense  research  activity  in  modelling  acoustic  propagation.  In 
particular,  in  the  language  of  the  sonar  engineer,  the  TL,  or  transmission 
loss,  of  this  energy  is  of  major  importance  for  i)  predicting  the  character  of 
the  sound  field  at  a  receiver  in  future  experiments,  ii)  for  comparing  crustal 
loss  with  the  better  known  TL  of  paths  remaining  primarily  in  the  water  layer, 
and  iii)  expanding  the  role  of  arrival  amplitudes  in  inversion  theory.  Just  as 
there  may  be  a  number  of  possible  paths  in  the  sea  between  a  source  and 
receiver,  each  with  a  different  loss  characteristic,  trajectories  in  the  crust 
are  variegated  and  exhibit  different  TL  behaviors.  It  is  important  to  be  able 
to  differentiate  the  energy  partitioned  among  the  different  paths,  and  to 
determine  which  paths  are  most  important. 

Resolving  the  locus  of  a  particular  acoustic  path  is  intimately  tied  to 
the  problem  of  determining  the  velocity  structure  of  a  medium.  To  the  limits 
of  the  geometrical  optics  approximation  of  acoustic  behaviour,  sometimes 
sorely  pressed  at  low  frequencies,  a  completely  detailed  knowledge  of  sound 
speed  variations,  both  laterally  and  with  depth,  plus  known  source 
characteristics  and  attenuation  losses  in  the  medium,  enables  one  in  principle 
to  predict  signals  observed  at  a  receiver.  For  an  ocean  acoustician,  the 
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requirement  of  environmental  knowledge  of  the  sound  speed  profiles,  both  in 
water  and  crust,  needed  to  predict  the  amplitude  and  timing  of  data,  is 
clearly  very  burdensome.  In  the  past  twenty  five  years,  however,  models  of  the 
oceanic  crust  have  been  formulated  which  are  statistically  consistent  over 
much  of  the  oceans.  These  models  divide  the  crust  into  three  or  more 
norizontal  layers  with  certain  average  thicknesses  and  velocities  (Raitt, 
1963).  At  least  within  the  confines  of  these  models,  if  a  typical  transmission 
loss  were  known  for  each  of  these  layers,  an  acoustician  can  make  predictions 
of  the  expected  strength  and  timing  of  crustal  arrivals  at  other  stations. 

Most  of  this  environmental  information  has  been  obtained  from  refraction 
and/or  wide  angle  reflection  data,  usually  via  travel  time  analysis.  Little 
has  been  done  in  developing  models  accounting  for  amplitude  dependence. 

Arrays  for  Refraction  Experiments 

The  standard  technique  in  ocean  refraction  experiments  has  basically 
involved  one  ship  and  one  or  more  receivers  (sonobuoy  or  OBS),  with  increasing 
range  between  shots.  With  a  dense  shot  spacing  and  large  enough  total  range, 
the  use  of  event  arrival  times,  especially  the  first  ones,  for  the  most 
prominent  features  in  the  data  has  been  sufficient  for  obtaining  a  reasonably 
good  understanding  of  the  velocity  structure  of  the  crust.  The  crustal  model 
referred  to  above  was  developed  from  averaging  experiments  of  this  type  from 
many  diverse  areas.  Lately,  a  multichannel  hydrophone  array  has  replaced  the 
single  receiver  in  some  experiments,  with  the  array  sometimes  being  towed  by  a 
second  ship  (Stoffa,  Buhl,  1979).  In  the  latter  technique,  termed  an  expanding 
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spread  profile  (ESP),  the  two  ships  start  at  a  common  point  and  steam  in 
opposite  directions.  In  this  way,  a  common  depth  point  is  shared  by  all  shots, 
with  the  use  of  a  Raydist  apparatus,  accurate  range  information,  which  is  a 
sensitive  parameter  in  the  inversion  methods,  is  also  available.  As  with  all 
arrays,  the  SNR  for  the  detection  of  coherent  energy  can  be  improved  with 
appropriate  processing.  Moreover,  estimates  of  the  received  energy  for 
different  horizontal  phase  velocities  can  be  made  which,  under  the  condition 
of  horizontal  crustal  layering,  provides  us  with  crustal  velocity  estimates 
using  just  one  shot.  However,  for  a  single  offset,  complete  information 
concerning  crustal  structure  is  not  be  obtained  since  the  SNR  for  certain 
events  is  range  dependent. 

Since  receiver  arrays  have  the  ability  of  generating  phase  velocity 
information  on  a  shot  by  shot  basis,  the  process  of  traveltime  analysis  used 
in  inversion  studies  can  be  somewhat  automated.  The  original  procedures  of 
generating  a  travel  time  versus  range  plot  for  a  sequence  of  densely  spaced 
snots  and  visually  picking  arrivals  can  be  improved  by  using  an  array  velocity 
analysis  technique  that  can  assign  velocities  to  arrivals  in  each  shot  trace. 
An  expanded  use  of  data  received  from  one  shot  would  minimize  interpretation 
errors  caused  by  uncertainties  in  range  and  source  level  variations.  Clearly, 
once  a  composite  of  a  number  of  shot  traces  is  developed  with  estimated  phase 
velocities  along  the  trace  for  each  shot,  the  problem  of  selecting  different 
arrival  times  for  a  particular  velocity  is  eased,  and  the  intercept  times  can 
be  found  for  use  in  traveltime  inversion  techniques,  eg.  the  tau-p  method 
(btoffa,  Diebold,  &  Buhl,  1981).  Array  processing  techniques  are  also 


important  in  discriminating  distinct  phenomena  that  occur  In  the  multipath 
reverberation  one  encounters  after  the  first  refracted  arrival,  and  effects  of 
local  inhomogeneities  such  as  bathymetric  variations  in  exploration  and/or 
oceanographic  experiments. 


Velocity  Analysis 

A  conventional  way  of  doing  array  velocity  analysis  employs  a  statistic 
that  estimates  the  amount  of  trace  to  trace  coherence  across  the  array,  for  a 
given  assumed  phase  velocity.  All  realistic  velocities  are  scanned,  and  the 
normalized  statistic,  a  "semblance  coefficient",  indicates  the  relative  amount 
of  energy  in  the  data,  at  each  velocity  (Sereda  and  Hajnal,  1976).  Another 
method,  used  throughout  in  what  follows  here,  employs  a  data  adaptive  spectral 
estimator.  Several  data  adaptive  techniques  were  originally  developed  in 
various  areas,  particularly  large  aperture  teleseismic  arrays  and  sonars.  The 
Maximum  Liklihood  Method  (Capon,  1969;  Edelblute,  1967;  Lacoss,  1971)  was  used 
at  Woods  Hole  originally  in  the  processing  of  reflection  data 
(Leverette, 1977) ,  and  eventually  extended  to  seismic  refraction  work 
(Baggeroer  and  Falconer, 1981).  The  technique  conceptually  designs  a  beamformer 
based  on  the  input  data  (hence  data  adaptive).  This  beamformer  minimizes 
output  power  with  the  constraint  that  energy  from  a  specific  direction  is 
passed  undistorted.  We  shall  see  that  the  structure  of  this  beamformer  can  be 
used  to  define  an  algorithm  that  estimates  what  is  known  as  the 
frequency-wavenumber  function  of  the  acoustic  field  for  a  certain  spatial 
frequency  associated  with  a  specific  direction.  Insofar  as  the  directions  of 
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the  arrivals  at  the  array  are  related  to  the  crustal  sound  speed  of  the  paths 
the  energy  has  traversed,  estimated  directions  lead  to  estimated  velocities, 
in  the  horizontal  layering  situation,  this  relationship  is  quite  simple  and 
the  velocities  estimated  are  very  accurate,  especially  at  high  SNR. 

In  stochastic  process  theory,  the  power  spectral  density  function  is  a 
measure  of  the  partitioning  of  energy  in  a  process  with  respect  to  frequency. 
The  corresponding  function  for  a  wide  sense  stationary  random  process  in  space 
and  time  is  the  frequency-wavenumber  function.  It  is  a  measure  of  the  mean 

:  square  power  per  unit  bandwidth  in  temporal  frequency  arriving  from  a  unit 

i 

l 

:  steradian  in  spatial  frequency  or  wavenumber,  which  is  uniquely  related  to 

horizontal  phase  velocity.  The  estimated  function  indicates  the  amount  of 
energy  that  has  arrived  at  the  array  via  a  particular  path. 

The  acoustic  field  generated  by  an  explosion,  however,  cannot  be 
modelled  as  a  stationary  process.  With  the  transient  nature  of  the  field,  only 
a  small  part  of  the  data  is  used.  This  "windowed"  data  must  then  be  treated  as 
if  it  were  part  of  an  ongoing,  time  invariant  process.  The  power  estimated  in 
tne  hypothetical  process  is  an  indication  of  the  actual  energy,  needed  for 
true  amplitude  measurement,  in  the  windowed  data  segments  that  were  employed. 
Tne  concept  of  windowing  data  to  track  nonstationary  phenomena  is  extensively 
used  in  signal  processing,  particularly  speech  analysis.  This  technique  is 
often  referred  to  as  "short  time,  spectral  estimation". 

The  MLM  estimate  is  known  to  be  biased  (Capon,  1969).  An  analytic 
expression  for  this  bias  has  yet  to  be  developed  for  all  possible  situations, 
however.  We  an  empirical  technique  that  can  be  used  to  evaluate  the  bias  for 


the  particular  data  set  and  array  configuration  discussed  below.  Given  a 
accurate  estimate  of  the  frequency/wavenumber  function  and  the  energy  spectrum 
of  a  source,  the  transmission  loss  for  a  certain  ray  path  can  be  determined. 
The  Kose  Experiment 

The  MLM  algorithm  and  our  transmission  loss  calculation  procedure  will 
be  applied  to  a  data  set  obtained  from  a  large  scale  acoustic/seismic  program 
(RUSE)  conducted  off  the  western  coast  of  Mexico  in  January  1979,  near  the 
East  Pacific  Rise.  Together  with  seventy  ocean  bottom  seismometers,  a  vertical 
(MARS)  and  a  horizontal  (ESP)  array  were  used  to  receive  acoustic  energy 
generated  by  a  series  of  explosions.  The  horizontal  array  was  towed  so  that 
data  was  received  in  the  ESP  format  described  above.  The  vertical  array  was 
stationary.  The  use  of  these  two  types  of  array  deployments,  and  of  the  bottom 
receivers,  resulted  in  one  experiment  employing  most  of  the  techniques 
currently  used  in  seismic  refraction  work. 

Insofar  as  the  experiment  occurred  near  an  active  plate  boundary,  the 
structural  makeup  of  the  crust  was  not  "typical",  and  difficulties  were 
experienced  in  relating  the  velocity  estimates  obtained  from  single  shots  to 
the  simple  layered  models  discussed  above.  As  we  shall  see,  the  complex 
seafloor  topography  also  limited  the  accuracy  of  our  calculated  velocities, 
however,  interesting  and  useful  results  were  obtained  and  estimates  of  crustal 


energy  partitioning  shall  be  presented. 
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Uverview 

In  Chapter  II,  a  summary  of  the  standard  theories  of  seismic  refraction 
is  given.  The  emphasis  is  on  current  ideas  concerning  the  strength  of 
refracted  waves.  Next  we  discuss  the  data  set  and  describe  the  different 
experiments  conducted  in  the  ROSE  project.  Chapter  IV  deals  with  the  velocity 
spectral  algorithm  and  the  method  used  to  determine  bias  corrections.  Chapter 
V  presents  some  results  of  the  computations  done  on  the  data  with  respect  to 
velocity  estimation.  Next,  we  describe  the  compensations  that  were  necessary 
to  make  the  measurements  obtained  from  the  algorithm  correspond  to 
transmission  loss  estimates  in  physical  units.  Source  levels,  biases,  surface 
effects,  group  beampatterns,  sensitivities,  and  analog  to  digital  conversion 
factors  must  all  be  included  to  arrive  at  estimates  of  path  losses.  Finally,  a 
summary  of  transmission  loss  estimates  from  this  data  set  is  presented. 
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CHAPTER  II 
SEISMIC  REFRACTION 

In  this  chapter,  the  fundamental  concepts  underlying  the  refraction 
experiment  are  presented.  In  particular,  we  concentrate  on  factors  influencing 
the  travel  time  and  amplitude  of  arrivals.  The  material  discussed  is  mainly  a 
review  for  the  geophysicist,  but  may  not  be  as  familiar  for  the  ocean 
acoustician. 

We  begin  with  the  free  space  solution  of  the  wave  equation  in  a 
homogeneous,  isotropic  elastic  solid.  We  then  discuss  acoustic  propagation  in 
a  simple  layered  medium,  with  one  interface  separating  two  isovelocity  half 
spaces.  Using  a  high  frequency,  ray  theory  analysis,  the  concept  of  a 
critically  refracted  interface  wave  is  presented.  We  show  that  this  analysis, 
based  on  the  "geometrical  optics"  model  of  sound  propagation,  does  not  explain 
empirical  observations  of  remotely  sensed  acoustic  events,  and  turn  to  a  "wave 
theory"  analysis  in  which  the  concept  of  "head  waves"  is  introduced.  Travel 
times  in  layered  media  are  accurately  predicted  by  head  wave  theory.  A  second 
interface,  representing  the  sea  surface,  is  added  to  the  model  and  we  define 
specific  events  observed  in  the  ROSE  data  which  can  be  represented  in  terms  of 
head  waves  and  surface  reflections  from  this  model.  Since  the  ocean  crust  is 
not  an  isovelocity  layer,  the  model  is  then  extended  to  include  multiple 
interfaces  Delow  the  seabed.  Events  received  at  different  horizontal  offsets, 
based  on  the  multiple  layer  model  and  head  wave  theory,  are  presented  in  the 
form  of  a  theoretical  travel  time/offset  (T-X)  diagram.  Because  of  the  absence 
of  events  that  correspond  to  expected  interlayer  reflections  with  this  model 
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in  most  refraction  data,  the  model  of  the  crust  is  finally  generalized  as  a 
region  with  a  continuous  velocity  gradient.  The  current  perspective  of  oceanic 
crust  is  based  on  this  last  model,  which  provides  better  agreement  with 
observed  arrival  amplitude  behavior.  We  show,  however,  that  some  layers  or 
interfaces  of  the  classical  layered  model  of  the  crust  mentioned  in  Chapter  I 
have  counterparts  as  regions  with  very  small  or  very  large  velocity  gradients 
respectively  in  the  continuous  model.  Finally,  since  we  are  concerned  with 
energy  partitioning  in  the  crust,  current  theories  of  head  wave  (in  layered 
media)  and  ray  (in  continuous  models)  amplitude  behavior  with  respect  to  range 
are  presented. 


Free  bpace  Propagation 

Let  ar  and  be  defined  as  the  Fourier  transforms  of  the 
dilational  and  rotational  displacement  potentials  in  an  elastic  solid.  Under 
the  conditions  of  homogeneity  and  isotropy,  the  Helmholtz  equations  in  free 
space  for  these  quantities  are  (Grant  &  West,  1965): 

V*  +*  Jk?  3l  — O  (2-la) 


(2-lb) 


wnere: 


"/vT^F-  <2-2» 

Jk0  -  ^  (2-2b) 

\  and  are  Lame's  constants,  to  is  the  radian  frequency,  and  ^  is  the 
density  of  the  solid.  In  a  homogeneous,  free  space,  two  dimensional  geometry, 
a  solution  found  by  separating  variables,  is  given  by: 

%  (*,»  cS)-  A  (oJ)  (JU- 


'*sa* 
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where  A  +  /*t  =  3-  ,  '£=  k^  /2ir,  and  A(  to),  corresponding  to  temporal 
benaviour,  is  an  arbitrary  function.  This  solution  represents  a  compressional 
(P)  plane  wave  traveling  with  velocity  in  a  direction  with  cosines 
{A  .0  »/*,)•  The  "wavenumbers"  kM  and  represent  spatial  frequency  in 
radians  and  cycles  per  unit  distance,  respectively.  The  solution  for  •xe  is 
the  same,  except  that  the  phase  speed  is  ^  ,  and  the  displacements  are 
orthogonal  to  the  direction  of  propagation,  representing  "shear"  (S)  waves. 

Medium  with  one  interface 

We  now  turn  from  the  free  space  model,  and  consider  a  medium  with  one 
horizontal  plane  boundary  separating  elastic  half-spaces  with  velocities 
andotj,  ^  as  in  Fig.  2-1  (Telford  et  al,  1976).  An  incident  compressional 
plane  wave  with  amplitude  h0  imposes  the  boundary  condition  that  apparent  wave 
numbers  in  a  direction  parallel  to  the  interface  are  constant.  This  leads  to 
Snell's  law: 

■^2—'  -  •^***■^1  X,  a 

'  Pa 

where  is  the  angle  both  of  incidence  and  of  P-wave  reflection,  (&*,At) 
are  the  angles  for  P  and  S  plane  wavefronts  that  are  "refracted"  into  the 
second  layer.  At  is  the  angle  of  reflection  for  an  S  wave  in  the  upper  layer 
and  the  constant  p  is  termed  the  ray  parameter.  If  the  sound  velocity  in  the 
second  layer  is  greater  than  ,  we  see  that  there  is  a  critical  incidence 
angle,  6^  ,  when  sin  6^  =  1.  At  incident  angles  ,  a  compressional  plane 
wave  solution  exists  that  travels  parallel  to  the  boundary  as  an  interface 


Geometry  of  reflected  and  refracted  waves  at  one  interface. 


Fig. 2-2 

Amplitudes  of  reflected  and  refracted  compressional  (P)  and  shear  (SV) 
versus  angle  of  incidence  at  one  interface  for  /  3  .33  and  / 

(from  Grant  &  West,  1965) 
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wave.  This  solution  Is  the  basis  of  all  simple  refraction  theories  and 
formulae.  The  "critically  refracted"  wave  travels  with  the  higher  speed.oC^, 
so  that  at  large  horizontal  ranges,  it  should  be  the  earliest  arrival. 

From  this  ray  theory  or  geometrical  optics  viewpoint,  however,  the 
interface  wave  will  not  appear  in  the  upper  layer,  and  its  predicted  amplitude 
is  zero.  The  latter  fact  is  seen  by  applying  six  boundary  conditions  of 
continuity  of  stress  and  displacement  at  the  interface  to  eqs.  2-1,  whereby 
the  Knott  equations  (Telford, 1976)  in  terms  of  the  potential  function 
amplitudes,  or  the  Zoeppritz  formulae  for  the  displacement  amplitudes  (Grant  & 
inest,  1965)  are  derived.  An  example,  calculated  from  these  equations,  is  shown 
in  fig.  2-2  (from  Grant  &  V/est)  in  which  ratios  of  incident  amplitude  to  the 
refracted  P  and  S  amplitudes  in  the  lower  layer  and  to  the  reflected  P 
amplitude  in  the  upper  layer  are  shown  versus  angle  of  incidence  for  a 
fluid-solid  boundary.  In  these,  ^*1/3  and  The  critical  angle 

for  the  compressional  (P)  and  shear  (SV)  waves  are  thus  sin  (1/3)  =  19.5°  , 

and  sin  (.6)  =  37°  ,  respectively.  In  the  ROSE  experiment,  typical  critical 

«  0 

angles  for  P  waves  were  in  the  10  to  15  range.  Note  that,  in  these  figures, 
all  energy  in  the  upper  layer  is  either  incident  or  is  reflected  from  the 
interface  at  angles  other  than  critical,  while  amplitudes  associated  with 
interface  waves  at  the  critical  angles  is  zero.  It  is  observed,  however,  that 
significant  energy  with  travel  times  much  as  one  would  calculate  for  an 
interface  wave  with  speed  refracting  energy  into  the  water  at  the  critical 
angle,  does  appear  in  refraction  experiments. 
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nead  waves 

The  most  widespread  theory  to  explain  this  is  based  on  Huygen's 
principle  using  curved  instead  of  ideal  planar  wavefronts.  It  predicts  the 
existence  of  "head  waves"  as  shown  in  Fig.  2-3  (from  Cerveny  and  Ravindra, 
1971).  In  2-3-a,  a  spherical  wavefront  originating  at  Me  impinges  upon  the 
interface  for  time  t *> h/oct  .  At  the  boundary  it  sets  up  a  disturbance  along  OP 
and  creates  Huygen  wavelets  (Fig  2-3-b)  which  produce  the  reflected  and 
refracted  wavefronts  where  constructive  interference  occurs.  The  speed  along 
the  interface  of  P,  the  contact  point  with  the  incident  wavefront,  is 

the  angle  from  PM.  to  the  horizontal  axis.  Beyond  a 

*  ‘/a. 

critical  horizontal  distance,  Xc  =  h  /  ((ot^/otj  )  -  1  )  ,  the  speed  of  this 

point  becomes  less  than  .  At  this  range,  the  angle  £(P)  has  increased  to 
the  critical  angle  .  We  now  get  the  situation  in  fig.  2-3-c.  The 

refracted  wave  in  layer  2  is  now  ahead  of  the  incident  or  reflected  fronts. 
Again  using  a  Huygen  construction,  M  1)  is  seen  to  be  a  locus  of  constructive 


interference,  and  for  constant  t*,  and  o*x,  is  a  straight  line  (in  2 
dimensions).  In  time  At,  the  disturbance  at  point  0*  will  move  both  to  Q 


along  the  boundary  at  speed  ,  and  to  point  M*  in  layer  1  at  speed  ot,  .  The 
angle  of  this  "head  wave"  is  seen  to  be:  sin"* =  6^  .  Wave  theory 
thus  predicts  that  the  Snell's  law  interface  wave  constantly  reflects  energy. 


in  the  form  of  a  head  wave,  back  into  the  uppermost  layer  at  the  critical 


angle.  The  apparent  horizontal  phase  velocity  of  the  head  wave  in  layer  1  will 


be: 


Fig. 2-3 

Construction  of  head  wave  using  Huygen  wavelets, 
(from  Cerveny  &  Ravindra,  1971) 


Thus,  if  the  direction  of  the  refracted  energy  or  the  horizontal  phase 
velocity  of  an  emerging  plane  wavefront  can  be  determined  within  the  upper 
medium,  the  sound  speed  in  layer  2  can  be  found  remotely  for  a  horizontally 
layered  medium.  This  is  the  basis  of  classical  inversion  theory  for  two  simpl 
layers  as  discussed  in  twing,  1963. 

Two  interface  model 

We  now  modify  the  proceeding  model  by  introducing  a  perfectly 
reflecting  interface  in  the  upper  half  space,  representing  the  sea  surface. 
Uue  to  surface  reflections,  many  arrivals  other  than  those  from  emerging 
interface  waves,  will  occur  at  a  receiver  with  this  model.  Referring  to  Fig. 
2-4,  together  with  the  critically  refracted  compressional  wave  labelled  IP,  a 
converted  shear  wave  (IS),  a  direct  wave,  and  a  series  of  water  layer 
reflections  (lw,  2W,  etc),  will  be  recorded  at  the  array.  Since  IP  refracts 
energy  continuously  back  into  the  water,  a  surface  receiver  may  encounter 
energy  which  travels  as  an  interface  wave  and  refracts  into  the  water.  Upon 
reflection  from  the  surface,  this  energy  reenters  the  seabed,  again  as  an 
interface  wave,  before  finally  refracting  into  the  water  and  being  detected. 
This  "multiple  refraction"  is  termed  2P  in  Fig.  2-4.  More  multiples  of  this 
type  (3P,  4P,  etc),  for  which  an  arrival  has  had  a  number  of  encounters  with 
the  surface,  can  be  observed  with  velocity  analysis  in  the  array  data 
presented  subsequently.  Often  it  is  found  that  the  amplitudes  of  2P  arrivals, 
and  even  those  of  higher  multiples,  are  stronger  than  IP.  Since  a  multiple 
refraction  arrival  can  be  the  sum  of  a  large  number  of  rays  each  travelling 


Fig. 2-4 

Possible  ray  paths  with  a  two-interface  model. 
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along  a  different  path,  this  is  possibly  due  to  the  constructive  interference, 
while  the  exact  acoustics  solution  still  remains  unsolved  for  these  multiples, 
they  are  important  from  the  ocean  acoustician's  perspective  because  of  their 
relatively  high  energy  levels. 

Multiple  layers 

The  upper  layer  in  the  two  interface  model  discussed  above  represents 
the  water.  For  the  region  below  the  seabed,  we  first  introduce  a  multiple 
layer  model  and  use  head  wave  theory  to  predict  events  received  at  a 
horizontal  offset  X  in  the  form  of  a  travel  time/offset  (T-X)  plot.  As 
discussed  in  Chapter  I,  the  original,  "classical"  model  of  the  oceanic  crust 
has  3  isovelocity  layers  above  the  mantle  interface. 

In  a  multiple  layered  case,  the  number  of  events  one  can  expect  is 
large,  especially  in  sedimentary  locales.  Figure  2-5  depicts  a  situation  with 
Im  interfaces.  In  constructing  a  time  versus  range  (T-X)  plot  for  this  example, 
and  concentrating  only  on  critical  refractions  of  first  arrivals,  we  see  that 
up  to  range  Xfit  ,  the  first  event  is  the  "IP"  from  the  layer  with  velocity  Vw. 
When  the  range  exceeds  X0 ,  the  IP  event  from  the  second  boundary  arrives 
earlier.  With  densely  spaced  sample  points  in  range,  the  locus  of  the  first 
arrival  traces  a  straight  line  in  the  T-X  plane  with  slope  1/V,  .  As  range  is 
increased  beyond  XCj  ,  the  interface  wave  from  the  Vz layer  will  eventually  be 
the  earliest  arrival.  This  pattern  continues  until,  at  the  largest  distance, 
tne  slope  of  the  first  arrival  line  will  be  1/VM  .  In  this  way,  for  a 
horizontal  layered  situation,  in  which  layer  velocities  increase  with  depth. 


Theoretical  T-X  diagram 
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tne  calculation  of  slopes  from  earliest  arrivals  on  a  T-X  plot  is  sufficient 
tor  obtaining  the  layer  velocities  of  interest. 

For  the  multiple  layer  case,  we  have  discussed  arrivals  due  to  head 
waves  only.  A  T-X  diagram  for  the  multi-layer  model  is  more  complex  than  this 
because  multiple  reflections  and  multiple  refractions  from  each  of  the 
interfaces  are  present.  These  appear  after  the  first  arrival.  Except  for  those 
involving  the  water  surface,  seafloor,  sedimentary  layered  sequences,  and  the 
mantle  interface,  interlayer  reflections  are  rarely  seen  (Ewing  &  Houtz,  1969) 
in  refraction  data,  however.  An  example  of  an  actual  T-X  plot  is  shown  in 
figure  2-6a  (Uetrick  &  Purdy,  1980)  for  an  experiment  conducted  near  the  Kane 
fracture  zone.  The  locus  of  events  that  can  be  attributed  to  layer  reflections 
are  limited  to  those  designated  by  PmP,  PmPPmP,  and  SmS,  for  the  mantle 
interface,  and  PnM  for  the  ocean  surface  boundary  (see  key  to  path 
nomenclature  in  fig.  2-6b).  In  the  ROSE  data  presented  subsequently,  which  is 
for  a  young  area  with  little  sediment  coverage,  the  only  clearly  identifiable 
reflections  we  find  involve  the  water  layer,  directly  (1W,  2 W),  or  indirectly 
(2P,etc.).  Since  strong  reflections  occur  at  areas  of  considerable  contrast  in 
elastic  properties,  e.g.  at  the  interfaces  in  the  model,  the  lack  of  reflected 
energy  argues  against  clearly  defined  layering  in  the  sub-basement.  Because  of 
tnis  experimental  evidence,  a  model  based  on  a  continuous  velocity/  depth 
relationship  in  the  crust  is  more  appropriate,  although  more  complex, 
especially  in  sedimentary  regions. 
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Continuous  velocity/  depth  profiles 

In  the  analysis  of  acoustic  propagation  in  regions  with  continuous 
velocity  gradients,  ray  theory  is  the  primary  tool  for  analysis.  Head  wave 
tneory  is  not  applicable.  Whereas  is  a  layered  model,  the  horizontal  phase 
velocity  of  an  emerging  plane  wavefront  is  equal  to  the  layer  velocity,  in  ray 
tneory,  it  is  equal  to  the  velocity  at  which  the  ray  turns  upward.  This 
velocity  is  also  the  reciprocal  of  the  ray  parameter,  i.e.  1/p.  Knowledge  of 
arrival  direction  remains  important  as  a  tool  for  remotely  obtaining 
information  about  velocity  structure. 

tty  conceptually  allowing  layer  thicknesses  in  the  multiple  layer  model 
to  approach  zero,  we  can  see  that  any  arbitrary  velocity/depth  relationship, 
provided  that  velocities  increase  with  depth,  may  be  determined  from  the 
slopes  produced  by  the  first  arrival  events  on  a  T-X  diagram.  This  is  the 
basis  of  the  Herglotz-Wiechert  integrals  (Aki,  1980),  used  for  the  inversion 
of  travel  time  versus  range  data.  The  density  in  range  with  which  data  is 
available  is  of  critical  importance  since  the  slope  of  the  first  arrival 
changes  continuously  with  offset. 

The  tau-p  method  of  travel  time  inversion  (Stoffa,  Diebold,  &  Buhl, 
1^81)  assumes  a  continuous  velocity  gradient  in  the  crust.  Using  first  arrival 
times,  it  produces  bounds  on  possible  velocity-depth  profiles  which  are 
compatible  with  the  data.  Many  of  these  resulting  profiles  are  similar  to  the 
example  in  figure  2-7  (Kennett,  1977).  Although  a  continuous  gradient  is 
assumed,  a  consistent  feature  of  these  velocity/depth  profiles  is  an  area  of 
small  velocity  change  immediately  above  a  sharper  gradient  at  the  mantle 
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interface.  This  region  corresponds  to  the  "classical"  layer  3  in  Raitt  (1963). 
Averaging  from  many  experiments,  this  region  has  a  mean  compressional  speed  of 
0.8  km/sec,  a  thickness  of  about  5  km,  and  begins  at  a  mean  depth  of  2  km 
oelow  basement.  We  shall  see  that  velocities  in  this  band  were  the  most 
prevalent  at  the  arrays  in  the  Rose  experiment.  Rays  that  travel  within  this 
"layer"  emerge  as  first  arrivals  at  offsets  of  approximately  10  to  30  km  in 
areas  where  ocean  depth  is  on  the  order  of  3  km.  Beyond  about  30  km,  mantle 
reflections  and  mantle  interface  waves  appear  as  the  earliest  events. 

The  model  of  the  oceanic  crust  we  have  been  employing  has  evolved  from 
a  simple  one-interface  case  to  a  continuous  velocity  gradient  representation. 
Although  the  latter  is  the  most  general,  we  point  out  that  at  least  the  mantle 
interface  and  the  sea  surface  can  be  effectively  considered  from  the  simpler, 
layered  model.  Reflections  from  these  interfaces  are  routinely  observed  in 
refraction  work  and  the  concept  of  head  waves  predicts  travel  times  along  the 
mantle  interface  accurately.  In  many  instances,  "layer  3"  can  also  be  treated 
as  a  homogeneous  isovelocity  layer. 

Amplitude  considerations 

Since  the  main  focus  of  this  paper  centers  on  energy  partitioning  in 
the  crust,  we  now  look  at  some  theories  concerning  head  wave  amplitude 
behavior  with  range  (for  layered  models)  and  the  amplitude  behavior  of  rays, 
(for  the  continuous  velocity  gradient  case). 

The  behavior  of  head  wave  amplitude  with  range  is  a  controversial 
issue.  Cerveny  and  Ravindra  (1971)  use  a  first  order  ray  series  solution  in 
solving  the  equations  of  motion  for  a  single  interface  problem  in  contrast  to 
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the  above  "zeroeth"  order  plane  wave  or  geometric  optics  solution.  They  obtain 

-3/1  -‘/j. 

an  amplitude  distance  curve  for  a  "pure"  head  wave  that  behaves  as  L  X 
wnere  X  is  the  horizontal  offset,  and  L  ■  (X-Xc)  is  the  propagation  distance 
along  the  interface.  According  to  this  equation,  at  large  ranges,  amplitude 
decreases  as  l/X1  ("spherical  spreading"). 

Alternatively,  if  we  assume  a  velocity  distribution  that  varies 
arbitrarily  with  depth,  ray  theory  predicts  that  the  pressure  amplitudes  will 
behave  with  distance  as: 


Pc’-Ro^e  „ 

my  «■. 
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(from  Clay  &  Medwin  (1977)) 

Re  is  the  reference  range  where  sound  speed  is  c# .  6&is  the  initial  angle  of 
a  ray  bundle  of  width  and  amplitude  P®.  0  is  the  average  angle  of  the 

bundle,  with  vertical  height  h,  at  horizontal  range  X  where  average  sound 
speed  is  c  (see  fig.  2-8).  This  equation  is  valid  at  ranges  where  focusing  of 
the  ray  bundle  does  not  occur,  so  that  Q  closely  approximates  the  angle  of 
all  rays  in  the  bundle  at  X.  For  a  source  at  the  surface  of  an  isovelocity 
layer  lying  above  a  half-space  with  a  linear  sound  speed  gradient  with 
slope  =b,  the  equation  for  the  mean  square  pressure  at  the  surface  at 
offset  X  reduces  to: 


pi  _  Jr  Po* 
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Rays  will  behave  with  an  x"1  amplitude  dependence  for  a  linear  gradient.  This 
type  of  geometrical  behavior  is  termed  "cylindrical  spreading"  and  is  also 


discussed  in  Kennett  ( 1977)  specifically  with  regard  to  crustal  acoustics.  We 
subsequently  see  that  TL  calculations  done  on  the  ROSE  data  yield  results 
wnicn  suggest  an  amplitude  attenuation  that  increases  somewhat  faster  than  the 
X  dependence.  This  is  probably  due  to  geometrical  losses  of  the  types 
discussed  discussed  above,  coupled  with  absorption  losses  in  the  crust. 

We  have  mentioned  that  at  an  approximate  offset  of  30  km,  energy  which 
nas  interacted  with  the  mantle  overtakes  "layer  3"  energy  as  the  earliest 
event.  Ray  bundles  with  different  parameters,  p,  appear  at  the  same  offset  and 
interact  to  produce  a  "focusing"  effect,  so  the  measured  transmission  loss  at 
these  ranges  will  be  low.  Although  TL  may  obey  a  cylindrical  or  spherical 
relationship  with  respect  to  range  from  an  overall  point  of  view,  fine  scale 
behavior  can  depart  from  the  general  trend  at  certain  offsets.  The  effects  on 
Tl  at  the  30  km  offset  will  be  shown  in  Chapter  5. 

Although  the  ROSE  data  to  be  presented  was  not  sampled  with  sufficient 
density  for  a  detailed  crustal  velocity  analysis,  events  from  "layer  3"  and 
the  mantle  interface  are  observed  with  the  use  of  the  velocity  analysis 
routine.  With  the  MLM  algorithm,  events  occuring  after  the  first  arrival  were 
also  identifiable.  We  report  on  these  results  following  a  description  of  the 
ROSE  experiment  and  a  discussion  of  the  analysis  algorithm. 
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CHAPTER  III 
THE  ROSE  DATA  SET 


The  Rivera  Ocean  Seismic  Experiment  (ROSE)  was  a  large 
seismic/acoustic  program  involving  ten  oceanographic  institutions  and  Navy 
Laboratories.  Originally  planned  to  be  sited  near  the  Rivera  Fracture  Zone,  it 
was  relocated  to  an  area  north  of  the  Clipperton  Fracture  Zone  because  of 
difficulty  in  obtaining  permission  to  operate  in  Mexican  territorial  waters. 
The  experiment  took  place  in  the  first  two  months  of  1979.  Figure  3-1  shows 
tne  general  area  of  the  experiment,  and  Fig.  3-2  maps  the  locations  of  some  of 
tne  instruments  deployed  with  respect  to  the  East  Pacific  Rise  central 
anomaly.  Seventy  ocean  bottom  seismometers,  a  12  channel  vertical  array 
(MAbS),  and  a  24  cnannel  towed  array  were  used  in  conjunction  with  explosive 
sources  ranging  in  weight  from  .1  to  1000  kg.  Five  research  vessels  were 
involved  in  the  project.  The  experiment  was  designed  to  study  the  following 
problems: 

1)  structure  and  evolution  of  young  oceanic  crust, 

2)  structure  and  dynamics  of  the  East  Pacific  Rise, 

3)  structure  and  dynamics  of  the  Orozco  fracture  zone, 

4)  long  and  short  range  propagation  of  low  frequency  acoustic  energy, 

b)  partitioning  of  energy  transmission  between  the  ocean  volume  and 

tne  crust/ lithosphere. 

To  investigate  these  problems,  the  experiment  was  divided  into  two 


phases.  Phase  I  consisted  of  an  active  program  of  shooting  explosives  to 
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the  bottom  instruments  and  to  the  acoustic  arrays.  Phase  II  was  primarily  a 
passive  earthquake  listening  experiment  with  some  calibration  shots.  At  Woods 
Hole,  Ur.  b.  Michael  Purdy  has  been  involved  with  processing  and  interpreting 
our  UBS  data,  while  Kenneth  Prada,  Thomas  O'Brien,  and  David  Sever  of  the 
Signal  Processing  Group  have  received  and  worked  on  data  from  the  ESP 
experiments  that  were  recorded  on  both  the  vertical  (MABS)  and  horizontal 
(ESP)  arrays.  Figure  3-3  shows  the  tracks  of  the  ESP  lines  that  were  shot,  and 
the  location  of  the  MABS  array.  Each  ESP  line  was  a  two-ship  experiment  with 
the  shooting  and  receiving  vessels  steaming  away  from  each  other  and  from  a 
common  midpoint.  Figure  3-4  is  a  schematic  of  the  shooting  schedule  used  in 
each  of  the  lines.  We  have  been  primarily  concerned  with  lines  2L  and  2S,  in 
which  the  midpoint  of  the  ship  tracks  was  in  the  near  vicinity  of  the  MABS 
array. 

The  Vertical  Array  (MABS) 

The  configuration  of  the  12  channel  vertical  array  is  shown  in  Fig. 
3-5.  The  data  tapes  were  recorded  in  analog  form  and  a  few  transcription 
difficulties  were  encountered.  The  MABS  required  high  amplifier  gains  for 
faithfully  recording  weak  refracted  arrivals.  The  subsequent  1W  arrivals  were 
then  of  sufficient  strength  to  saturate  the  analog  recorders,  with  their 
limited  dynamic  range.  Our  analysis  of  the  MABS  data  therefore  was  confined  to 
the  ("refracted")  arrivals  that  appear  before  overload.  Also,  hydrophones  5,7, 
and  11  only  worked  intermittently,  while  the  deepest  sensor,  12,  did  not 
function  at  all.  Digitization  of  the  data  was  necessary  for  use  in  the 
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velocity  analysis  programs.  Using  an  analog  tape  recorder,  the  MASS  tapes 
tapes  were  played  back  to  the  WHOI  digital  seismic  acquisition  system  (Prada 
et  al,  1 y 74) .  Due  to  problems  involving  the  synchronization  of  the  analog  tape 
speed,  the  recorded  time  codes  on  the  analog  tapes,  and  the  resulting  time 
information  on  the  digitized  tape  headers,  preliminary  velocity  analysis  for 
MABb  data  had  to  be  based  on  relative  times  calculated  from  direct  water 
arrivals  and  the  possibly  inaccurate  vessel  positions  of  the  ship  logs.  This 
navigation  information  was  based  mainly  on  a  SATNAV  system  which,  because  of 
low  latitudes,  was  only  updated  on  the  order  of  once  every  90  minutes  in  an 
area  of  strong  ocean  currents.  Estimated  position  errors  were  greater  than  1 
nautical  mile.  The  output  of  the  acquisition  system  was  in  WHOI's  12  channel 
CANBARX  format,  and  this  was  translated  to  ROSE  format.  The  CANBARX  format 
with  the  8  good  data  channels  was  used  in  the  MLM  processing  of  Line  2S  MABS 
data. 

The  Horizontal  Array  (ESP) 

Better  results  were  experienced  with  the  towed  array.  In  line  2S,  the 
shooting  ship  was  the  University  of  Hawaii's  R/V  Kana  Keoki.  It  dropped  5  and 
25  lb.  charges  as  it  steamed  away  from  the  area  of  the  vertical  array. 
Lamont-Uoherty's  R/V  Conrad  towed  the  seismic  streamer  array.  Each  active 
section  in  this  array  was  100  meters  in  length,  consisting  of  two  50  meter 
hydrophone  groups  connected  in  parallel,  with  no  taper.  The  total  length  of 
the  array  was  2400  meters  and  all  24  channels  of  the  data  were  good.  The  tapes 
were  received  in  Lamont's  digitized  field  format  and  were  translated  to  the 
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CANBrtrtX  format  for  processing  with  the  MLM  algorithm.  At  the  time  line  2S  was 
received,  existing  computer  memory  size  limited  the  processing  to  the  use  of 
only  twelve  channels  of  data  (see  Fig.  3-6).  By  the  time  line  2L  data  arrived, 
tne  incorporation  of  an  FPS  AP120B  array  processor  and  the  expansion  of  the 
system  permitted  working  with  all  24  channels  in  a  faster  version  of  the 
analysis  routine.  Data  was  translated  from  Lamont  to  the  SEGY  format,  which 
nas  Decome  the  standard  at  WHOI.  Unfortunately,  there  were  some  problems  with 
the  ESP  data  as  well.  Surface  reflections  at  low  frequencies  attenuate 
arrivals  from  directions  parallel  to  the  streamer,  due  to  the  LLoyd  mirror 
effect,  so  direct  water  waves  cannot  be  seen  with  the  array  processing 

procedure;  however,  the  water  bounce  arrivals  (1W,2W)  can  be  analysed. 

The  time  information  given  in  the  digitized  ESP  data  does  not  include 
fractions  of  a  second,  so  that  an  error  of  one  second  is  possible  in  travel 
times  based  on  this  information.  Time  estimates  based  on  direct  water  arrivals 
were  therefore  more  accurate  for  both  the  MABS  and  the  ESP  data.  This  is  the 
procedure  usually  followed  in  travel  time  and  range  calculations,  and  is  the 
oasis  for  travel  time  in  some  T-X  plots  presented  subsequently.  For  each  ESP 
shot,  40  seconds  of  data  were  sent  to  WHOI.  With  this  data,  estimates  for  all 

of  line  2S  were  made.  Because  of  the  larger  extent  of  the  2L  line  (beyond  100 

km.),  however,  water  arrivals  beyond  the  range  of  78  km.  were  not  available. 
Travel  time  estimates  were  therefore  not  made  for  the  latter  part  of  2L. 

we  have  shown  that  the  number  of  channels  of  data  used  in  processing 
KOSE  data  varied  as  different  situations  evolved.  The  bias  in  the  estimated 
received  energy  calculated  with  the  MLM  programs  is  dependent  on  the  number  of 
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24  channel  horizontal  array  geometry  and 
location  of  sensors  used  in  12  channel  processing. 
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channels  of  data  that  are  used.  So  that  consistent  transmission  loss  estimates 
would  be  obtained  from  sets  of  ROSE  data  processed  with  different  array 
configurations,  programs  designed  for  studying  MLM  bias  were  run  for  model 
arrays  with  8,  12,  and  24  data  channels.  These  are  discussed  in  Chapter  4, 
following  a  description  of  the  MLM  algorithm. 


CHAPTER  IV 


FREQUENCY -WAVENUMBER  FUNCTION  ESTIMATION 
This  chapter  discusses  estimation  of  energy  partitioning  at  the  arrays 
in  the  RUSE  experiment.  Waveform  characteristics  of  in  situ  acoustic  sources 
and  paths  are  not  known  in  detail,  so  that  a  received  waveform  is  modelled 
prodabalistical ly.  We  define  a  set  of  basis  functions  from  the  theory  of 
space/tiine  random  processes,  emphasizing  the  "frequency/wavenumber"  function, 
P(f,^).  This  function  is  based  upon  a  mathematical  representation  of  a 

spatially  homogeneous,  temporally  stationary  random  process  as  a  superposition 

i  2.-*  (  jLt-v--  &) 

of  plane  waves, C  5  ,  with  temporal  frequency  f  and  wavenumber,  V*. 

we  snow  that  P(f,V[)  is  a  measure  of  the  partitioning  of  energy  with  f  and  vf, 
in  a  process. 

In  applying  random  process  theory  to  propagating  acoustic  waves,  a 
constraint,  the  "dispersion  relation",  is  imposed  upon  the  relationship  of  f 
and  'Z  for  plane  waves:  W\  =  '/X  =•£/<!  ,  where  c  is  the  sound  speed,  and  X 
the  wavelength.  The  unit  vector  V"  /  1^1  is  in  the  direction  of  propagation.  In 
Chapter  II,  we  discussed  the  fact  that  knowledge  of  the  spatial  distribution 
of  coherent,  propagating  energy,  having  interacted  with  the  crust,  can  lead  to 
crustal  velocity  estimates.  Given  a  known  source  level  and  a  valid  crustal 
model,  we  also  noted  that  knowledge  of  the  magnitude  of  energy  arriving  at  a 
certain  spatial  angle  at  a  receiver  can  be  used  to  determine  the  transmission 
loss,  TL,  of  the  crustal  path  corresponding  to  the  angle.  Specifically,  if 
there  is  an  arbitrary  velocity/depth  relation  in  the  crust,  the  vertical  angle 
of  a  ray,  ©  ,  is  related  to  the  horizontal  phase  velocity,  c^,  of  the 
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propagating  energy  at  the  receiver  and  to  the  sound  speed  in  the  crust,  c^,  at 
which  the  ray  turns  upward: 

■r  -  --  'K  -  'A*  l4-’> 

where  p  is  the  ray  parameter.  If  the  acoustic  field  in  an  experiment  is 
modelled  as  a  homogeneous,  random  process,  the  estimation  of  P(f,  </)  for  a 
range  of  orientations  of  the  vector,  ✓  /  |^|  ,  generates  information  about 
the  directional  character  of  energy  in  the  field.  This  is  used  for  crustal 
velocity  estimation.  Likewise,  the  magnitude  of  the  estimated  wavenumber 
function  is  the  basis  for  TL  estimates.  However,  since  impulsive  sources  are 
used  in  most  refraction  experiments,  the  received  data  is  not  stationary.  By 
employing  short  segments  of  data  that  are  treated  as  samples  of  a  hypothetical 
stationary  process,  the  concept  of  the  frequency/wavenumber  function  is  made 
applicable  to  the  ROSE  data. 

The  estimation  of  the  frequency/wavenumber  function  is  done  with  an 
array,  which  essentially  measures  the  apparent  phase  velocity  of  coherent 
waveforms  along  its  geometry.  The  estimated  phase  velocity  leads  to 
directional  information  as  seen  above.  From  a  mathematical  viewpoint,  the 
array  is  treated  as  a  deterministic  spatial  sampler  of  random  processes:  its 
geometry  and  temporal  frequency  response  are  parameters  that  can  be  adjusted 
to  produce  a  certain  "frequency  wave-vector  response"  function,  >$7(f,^). 

This  function  specifies  the  response  of  the  array  to  a  deterministic  plane 
wave,  the  design  of  this  response  being  termed  "beamforming".  For  estimation 
of  P(f,*£),  the  ideal  response  function  will  only  pass  energy  in  a  small 
spatial  angle  corresponding  to  a  narrow  wavenumber  or  "spatial  frequency" 
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oand,  v,  and  completely  reject  energy  from  other  areas.  The  commonly  used 
"conventional"  or  delay  and  sum  beamformer  is  first  discussed.  It  is  well 
known  that  this  conventional  processor  has  large  sidelobes  for  sparse  arrays 
with  a  small  number  of  sensors.  An  optimal  array  processor,  the  MLM 
beamformer,  is  then  introduced  which  minimises  this  sidelobe  effect.  The 
output  power  from  this  processor  is  the  basis  of  our  estimates  of  energy 
partitioning. 


Space/Time  Random  Processes 

Stochastic  processes  for  time  series  have  a  one-dimensional  index  set, 

\ 

e.g.  t,  in  x(t).  In  the  array  processing  problem,  the  index  set  may  increase 
to  four  dimensions  as  in  x(t,x,y,z),  or  x(t,r).  Most  of  the  concepts  involved 
relate  directly  back  to  time  series,  although  the  array  introduces  unique 
considerations.  In  particular,  for  a  zero  mean,  wide  sense  stationary  time 
series,  x(t),  that  is  an  input  to  a  pair  of  ideal  linear  bandpass  filters  with 
responses  Ht(f)  and  H^f),  as  shown  in  Figure  4-1,  the  mean  square  power  in 
the  output  process,  y  (t),  is  (from  Appendix  4-1): 

where  Q,  denotes  expectation,  W  is  the  bandwidth  of  the  filter  with  center 
frequency  f.  ,  and  S  (f )  is  the  power  spectral  density  function  of  the  input 
process: 


(4-2) 


(4-3> 

Since  n  can  be  made  arbitrarily  small,  S^f)  is  a  measure  of  the  mean  square 
power/unit  bandwidth.  The  cross  correlation  function  of  the  two  output 
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processes,  K^('r),  is: 

=  j  i4-4) 

-\£^({)H, ({)«'({)*£  {  =  g!!S^rffnt 

Since  y(  (t)  and  y  (t)  are  both  derived  from  x(t),  this  implies  that  disjoint 
Dands  of  a  stationary  times  series  are  statistically  uncorrelated.  The 
frequency  representation,  (S*(f)),  provides  a  powerful  and  possibly  simpler 
area  in  wnich  to  work  with  the  time  series. 

ke  now  examine  the  analog  of  S^Jf)  for  space/time  random  processes. 

The  following  functions  for  a  stationary  (in  time),  homogeneous  (in  space) 
process,  x(t,r),  are  defined: 

Space/Time  Correlation  function: 

(4-5) 

Spectral  Covariance  function: 

(4‘6> 

Note  that,  in  general,  the  latter  is  the  cross  spectral  density  function 
S  (f)  for  the  time  series  x  (t)  and  x,(t)  at  r  and  r-*r.  If  Ar  =  0  ,  it 

/*\XX  i  *  -  - 

is  just  the  spectral  density  function  of  the  time  series  at  r. 
Frequency/Wavenumber  function: 

where  y*  is  the  wavenumber  or  spatial  frequency.  The  estimation  of  the 
frequency/wavenumber  function  with  an  array  is  the  object  of  the  MLM 
algorithm.  Irfhen  based  on  refraction  data,  the  orientation  of  where  the 
estimate  of  this  function  is  large  indicates  the  direction  of  arrival  of  a 
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s i gn i t i cant  event  and  is  used  for  crustal  velocity  estimation.  This  function 
is  defined  for  stationary,  homogeneous  space/time  processes  only.  In  the 
analysis  of  the  ROSE  data,  we  model  the  acoustic  field  as  being  short  term 
stationary  and  homogeneous. 

The  frequency/wavenumber  function  is  the  parallel  of  the  spectral 
density  function  in  time  series.  To  see  this,  we  next  discuss  arrays  which 
correspond  to  the  linear  filters  (H,(f)  and  Hjf))  in  fig.  4-1. 

A  discrete  array  with  K  sensors  at  locations  remeasuring  x(t,rx)  is 
shown  in  Fig.  4-2.  Each  sensor  is  connected  to  a  linear  filter  g^  (t),  and  the 
outputs  of  all  filters  are  summed  to  produce  the  array  processor  output,  y(t): 


“  £  g  jL  C£)  *  OC  C 


(4-8). 


where  *  is  the  convolution  operation.  If  x(t)  is  a  simple  sinusoidal  plane 
wave  of  temporal  frequency  f,  and  wavevector  <S  :  x(t)=  ,  the 


output  y(t)  is: 


(4-9) 


De*7*** 

where  b^( f )  is  the  frequency  response  of  the  filter  at  r^.  The  output  of  the 
array  is  just  a  sinusoid  with  amplitude  and  phase  determined  by  ,  the 

frequency/wave  vector  response  function  of  the  array.  This  is  the  function 
that  determines  the  beam  pattern  of  the  array. 

Delay  and  Sum  Beamformer 

If  we  wish  to  look  in  the  direction  of  the  unit  vector  1*1  for 
plane  waves  with  speed  c,  a  reasonable  inpulse  response  for  the  linear  filters 
in  fig.  4-2  would  be: 

(4_10) 
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Unly  the  "correct"  plane  waves  arriving  at  each  r.  sum  in  phase,  hence  the 
name  "delay  and  sum  beamformer"  for  this  particular  processor.  The  dispersion 
relation  constrains  the  region  in  the  frequency  domain  where  real  plane  waves 
can  propagate,  and  hence  the  region  of  interest  for  the  frequency/wavevector 
response  function.  Incorporating  this  relation  in  eq.  4-9  for  the  case  of  the 


delay  and  sum  beamformer,  we  obtain  the  response  function: 

In  particular,  for  an  aperture  that  is  a  straight  line  segment  of  length  L,  in 

A 

the  direction  a  ,  the  response  function  becomes: 


where  fi^is  the  angle  between  the  target  and  the  array  directions,  and  &  is 
the  angle  between  the  actual  arriving  plane  wave  and  the  array.  Figures  4.3, 
4.4,  and  4.5  depict  beampatterns  for  cases  in  which  is  90*  ("broadside"), 
45  ,  and  U  ("endfire"),  respectively.  Notice  that  the  sidelobes  in  these 
figures  can  pass  an  appreciable  amount  of  energy  coming  from  directions  other 
than  that  desired,  and  that  the  main  lobe  can  be  quite  wide,  especially  at 
endfire. 


Array  Response  to  Random  Processes 


Having  discussed  the  response  of  an  array  to  a  plane  wave,  we  now  look 
at  its  behavior  in  a  homogeneous  and  stationary  random  space/time  field  with 
spectral  covariance  S^f  :*.&).  For  the  discrete  array  in  Fig.  4-2,  the 
spectral  density  function  of  the  output  time  series,  y(t)  is  (Appendix  4-2): 


(4-13) 
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Fig.  4.3  Beam  pattern  of  a  uniform  line  array, 
L/\=  3. 5.  @,  =  ±16.6°,  $£=±34.9°, 

Qy±  59.0°.  The  pattern  is  rotationally 
symmetric  about  the  z  axis, 

{ broadside  case ) 
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Fig.  4.4  Beam  pattern  of  a  line  array  with 
uniform  amplitude  and  linear 
phase -shift,  L/X  =  3.5  ,  &ta=  tt/4. 

=  83.1  °  and  24.9°,  02=  7.8°, 

^  =  -8,6°,  Q^=  - 25.8 ° ,  05=-46.2°. 
The  pattern  is  rotationally 
symmetric  about  the  z  axis 


Same  array  as  in  Fig. 4.4, 
but  with’  0tQ=  7r/2. 

0,  =  45.6° ,  02=25.4°,  ^3=8.2°, 
04=- 8.2°,  #5=  -25.4 0g  =  - 45.6°, 
07  =  -9O.O°  ( end  fire  case  ) 
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lf  we  form  a  vector  &-&(&■■■■  <GM  • 
trix  this  expression  can  be  written  as  the  qua 

S,C#-QTCSt(<)]&' 


and  a  KxK 


matrix 


this  expression  can  be  written  as  the  quadratic  form: 


(4-14) 


by  expressing  tne  covariance  function  as  the  Fourier  transform  of  the 


frequency/wavenumber  function,  this  can  also  be  expressed  as: 

(-{)  =  P(|*)  I  v  (4-,5» 

If  £){u  \[)  is  a  "pencil  beam"  response  function,  having  unit 
magnitude  over  narrow  spatial  and  temporal  frequency  bands,  V  and  W,  and  being 
zero  elsewhere,  the  array  is  the  counterpart  of  the  ideal  bandpass  filter  in 
fig.  4-1.  This  can  be  seen  by  evaluating  the  mean  square  output  power  of  the 
processor  with  this  response  function: 

^  =R^.(o)  =  1  df  SSS  <4-16a) 

WV  (4-,6b) 

This  equation  is  the  analog  of  eq.  4-2  for  time  series.  Conceptually,  V  and  W 


(4-16a) 
( 4- 16b) 


can  be  made  arbitrarily  small  so  that  P^f,*/)  represents  the  power  per  unit 
spatial  and  temporal  bandwidth  for  a  homogeneous,  stationary  process.  It  is 
the  spatial  analog  of  the  power  spectral  density  function.  For  a  homogeneous, 
stationary  process,  P(f,  ✓)  is  a  measure  of  the  energy  arriving  at  temporal 
frequency  f,  from  the  direction  represented  in  and  it  is  not  influenced  by 
energy  arriving  from  other  regions  in  the  frequency  domain. 

Optimal  response  function  for  arbitrary  noise 


Having  discussed  the  importance  of  the  function  P(f,  ✓),  we  now  turn 


to  the  problem  of  its  estimation.  From  eq.  4-16b,  we  see  that  the  power  at  the 
output  of  an  array  processor,  with  a  sufficiently  narrow  passband  with  unit 
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magnitude  in  temporal  and  spatial  frequency,  can  be  an  estimate  of  P( f ,  .  We 

require  a  j^7(f,\r)  that  is  unity  for  a  desired  "target"  ,  and  is  minimized 
in  all  other  frequency  regions.  In  Appendix  4-3,  an  optimum  processor  is 
determined  with  these  specifications  for  a  discrete  array  of  K  sensors  and  an 
aruitrary  stationary  process  with  covariance  S(f,  r^-rj)  =  S-.  (f).  In  matrix 
notation,  the  minimised  output  power  density,  subject  to  the  constraints 
mentioned  is: 


(4-18) 


C  (A ,  j- _ 

’  »  £Ttf4G*<pT6'({* I 

wnere  the  "steering  vector"  E(f,^t)  is: 

=  '/«  . 

S  (f)  is  the  estimator  for  P(f,<£).  For  a  Gaussian  space/time  process,  this  is 

9 

the  maximum- liklihood  (ML)  estimate.  The  expression  in  eq.  4-17  is  the  basis 
of  the  MLM  algorithm  used  in  the  ROSE  data  analysis.  Sidelobe  levels  and  null 
positions  in  the  optimal  >^7(f,^)  are  adjusted  so  that  noise  is  optimally 
attenuated  in  directions  where  interference  is  strong.  It  is  data  adaptive  in 

that  it  requires  knowledge  of  the  covariance  function  S^(f)  of  the  process 

o 

that  it  is  operating  on.  S  •• (f )  must  usually  be  estimated  from  the  data  for 

?  A 

each  implementation  of  the  estimator,  the  estimate  denoted  as  S--(f).  The  MLM 

<r 

estimator  for  the  frequency/  wavenumber  function  is  then: 


fL.tyc'i  * 


(4-19) 
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Me  now  discuss  two  examples  of  mathematical  process  models  for  which 
closed  form  expressions  for  the  covariance  can  be  written.  These  examples  are 
used  to  investigate  the  behavior  of  the  MLM  expression. 


uncorrelated  Sensor  Noise 

tven  in  situations  in  which  no  propagating  process  exists,  with 
factors  such  as  fluid  flow  noise,  thermal  noise,  etc,  a  sensor  output  is  never 
completely  deterministic.  A  random  term,  w(t,r),  often  considered  additive, 
will  be  present  at  each  sensor,  so  that  the  covariance  function  of  sensor 


outputs,  wlt,r),  is  often  modelled  as: 

(4-20> 

The  notation  SL  signifies  that  the  Kronecker  delta  is  defined  at  sensor 

t 

locations  only  and  that  noise  at  r^is  independent  of  noise  at  r^.  The 
process,  w(t,r),  is  not  homogeneous  in  space  and  cannot  be  described  by  a 
frequency-wavenumber  function.  From  eq.  4-13,  the  output  density,  $j[f),  for 
an  arbitrary  array  processor  in  this  noise  field  is  just: 

-  %  #, &{)  £.  ({)  % 

=  %  J{)\^  (4'2,b) 

in  Appendix  4-4,  the  optimal  response  function  for  this  particular  noise 


( 4-2 1  a) 
( 4-21 b) 


process  is  determined  to  be: 

h-22) 

This  processor  is  the  conventional  delay  and  sum  beamformer.  It  is  optimal  in 


this  case  of  spatially  uncorrelated  sensor  noise.  The  total  output  power 
density,  and  the  optimal  estimate  of  P(f,^),  is  then: 


•r- -<--*** 'Sifrwy; 


(4-23, 

In  estimating  P(f,V)  for  a  mathematical  model  in  which  its  value  for  all 
arguments  is  undefined,  this  estimator  returns  with  the  noise  power,  reduced 
by  a  factor  of  K,  the  “array  gain11. 

If  this  processor  is  used  as  an  estimator  of  P(f,vf_t)  in  a  situation 
in  which  the  process  covariance  is  arbitrary,  the  estimate  in  matrix  form  can 
be  written  as: 

L,  4^  •  <4-24» 

where  t(f,  is  the  steering  vector.  This  expression  is  often  known  as  the 
"conventional"  estimate  of  P  (f,  <c).  Recalling  the  beampatterns  associated 
with  the  delay  and  sum  response  function  in  figures  4-3  to  4-5,  however,  we 
see  that  significant  energy,  from  directions  other  than  \£,  is  passed  in  the 
sidelobes. 

Optimal  Estimation  of  Unidirectional  Process  in  Spatially  Uncorrelated  Noise 
The  second  example  of  the  application  of  the  MLM  estimator  is  a  case 
in  which  a  stationary  homogeneous  process,  x#(»),  propagating  from  one 
direction  only,  ^ /  t^l  ,  is  added  to  the  model  discussed  above: 

rL  C±,&)  =  O  -  -20-  (*:,■&)  (4-25> 

This  example  is  particularly  important  in  refraction  work,  since  an  arrival  is 

modelled  as  a  windowed  sample  of  this  type  of  random  process.  Assuming  in  this 

x  c  z. 

case  that  S^f.Ar)  ,  where  6  is  constant  with  respect  to  f  ("white 

noise"),  the  covariance  function  of  x(t)  is: 

-(4-> 


•  SW«: 
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wnere  P  (f)  is  the  power  spectrum  of  x0(*).  It  represents  the  power  in  the 
propagating  component  of  the  process.  In  this  model,  the  frequency/wavenumber 
function,  P(f,  *r)  is  equal  to  P  (f)  at  «/  =  >/"p ,  and  is  zero  elsewhere.  We 
now  employ  the  MLM  estimator  of  P(f,V).  For  a  discrete  array,  with  K  sensors, 
samples  of  the  total  covariance  function  at  sensor  locations  can  be  written  in 


matrix  form  as: 

K1  PJ{)  6Z I 

where  i  is  the  KxK  identity  matrix,  and  £(f ,  </p)  is  the  steering  vector.  This 
covariance  matrix  can  be  inverted  by  using  the  identity: 

(A*4i-iT)  *=  A'- *2:  A  1  (4-28) 

where  A  is  a  KxK  matrix,  u  and  v  are  KxM,  and  I,*,  is  the  MxM  identity  matrix. 

^  X  T 

We  identify  A  with  I,  E  (f,  ✓p)  with  u,  and  E  (f,  v^)  with  v  ,  so  that  M  =  1. 

If  we  substitute  this  into  the  expression  for  the  optimal  MLM  estimate  at 

T  * 

frequencies  f  and  and  recognize  that  E  E  =  1/K,  we  obtain: 

<  fcl  «  i.Mri-u‘1  -(4-29) 

where:  L  X  J 


.(4-29) 


Comparing  eq.4-3U  with  4-22,  we  notice  that  Q  represents  the  response  of  the 
conventional  array  processor  at  \T  when  the  target  is  \f  .  If  is  widely 
separated  from  ,  ^<atU,  and  P  (f,  •Q  is  just  6/k  ,  the  estimate  we 
found  in  a  sensor  noise  field  alone.  If  ^  =1,  and  we  get: 
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The  uncorrelated  noise  is  again  reduced  by  the  array  gain,  K.  With 
this  process  model,  the  MLM  estimate  in  the  direction  of  plane  wave 
propagation,  results  in  the  correct  value,  P^(f),  together  with  the  added 
noise  component  reduced  by  the  array  gain  K.  If  the  noise  power,  6*~  can  be 
estimated,  by  directing  the  array  where  no  coherent  energy  is  propagating,  the 
value  of  P  (f )  follows  immediately.  This  is  the  procedure  used  to  determine 
the  coherent  energy  arriving  from  a  particular  direction  in  a  refraction 
experiment,  after  data  is  suitably  windowed  so  that  this  model  is 
approximately  valid. 

implementation  of  the  MLM  Algorithm 

We  now  present  the  procedure  used  in  evaluating  eq.  4-19  for  the  ROSE 
data  set.  The  material  discussed  to  this  point  is  based  on  a  homogeneous, 
stationary  process  assumption.  Refraction  data  cannot  be  modelled  as 
stationary  since  it  changes  character  with  a  time  constant  determined  by  the 
source  signature.  The  essential  idea  behind  the  implementation  is  to  use  the 
stationary  concept  we  have  been  discussing  over  a  spatially  finite  and 
temporally  short  analysis  window.  For  an  estimate  at  frequencies  f  and  we 
model  the  data  as  samples  of  a  random  process  consisting  of  a  single 
propagating  plane  wave  at  wf^with  added  uncorrelated  noise,  as  in  the  example 
above.  This  is  usually  valid  in  refraction  work  if  the  data  segments  used  are 
short  enough  so  that  only  one  event,  corresponding  to  a  coherent  arrival  from 
one  specific  direction,  is  fully  represented  in  the  Windowed  data  from  all 
sensors.  In  the  implementation,  T  seconds  of  data  (typically  T  =  1  sec.  in  our 
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application)  are  used  from  each  channel  for  an  estimation  at  a  certain  travel 

time.  For  each  covariance  matrix  term,  S*.  (f),  and  target  direction,  K£j, 

w 

the  T  seconds  are  selected  for  each  pair  of  sensors  at  points  corresponding  to 
tne  times  the  hypothetical  plane  wave  would  appear  at  each  sensor.  The 
horizontal  phase  velocity  of  this  wave  is  related  to  V^by: 

~  ~  I  (4-32) 

where  a  is  the  vertical  angle  of  arrival  and  i i s  a  unit  vector  in  the 
vertical  direction.  The  data  within  the  shaded  "windows"  in  figs.  5.2  are 
examples  of  typical  segments  that  would  be  selected.  The  traces  are  received 
waveforms  from  12  of  the  24  sensors  of  the  horizontal  array  in  one  ROSE 
experiment.  In  fig  5.2a,  the  window  at  t=2  seconds  is  almost  horizontal,  i.e. 
with  little  "moveout".  Using  this  data,  the  estimation  procedure  models  the 
field  as  a  plane  unidirectional  process  with  a  relatively  high  horizontal 
phase  velocity  encountering  the  array  almost  at  broadside.  In  contrast,  in 
fig.  5-2b,  the  window  shown  indicates  a  relatively  larger  moveout,  so  that  the 
segments  are  considered  to  be  samples  in  time  and  space  of  a  plane  wave 
process  arriving  from  a  direction  closer  to  endfire. 

For  the  ROSE  data,  estimates  of  P(f,  V)  were  made  for  vertical  target 
angles  of  8*  to  90°,  and,  for  ESP  experiments,  at  the  azimuthal  angle  directly 
behind  the  receiving  ship.  In  a  horizontally  layered  crustal  model,  this 
corresponds  to  phase  and  layer  velocities  of  1.5  to  about  10  km/sec.  Once  T 
seconds  of  data,  with  the  proper  moveout,  is  selected  for  each  target,  the 
covariance  matrix  necessary  for  evaluating  eq.  4-19  is  formed.  Figure  4-6  is  a 
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block  diagram  of  the  procedure  used  for  one  matrix  term,  S- (f).  The  segments 
from  sensors  i  and  j  are  windowed  in  time  with  a  cosine  taper,  and  fast 
Fourier  transformed  (FFi  length:  N).  In  the  ROSE  implementation,  the  effective 
data  length  after  windowing  is  T  »  .5  sec.  After  taking  the  complex  conjugate 
of  the  coefficients  of  the  jth  sensor,  a  product  is  calculated  for  each 


coefficient. 


To  stabilize  the  covariance  estimate  at  a  frequency  f,  a  simple 
average  of  the  coefficient  products  over  a  band  of  width  W  Hz  centered  about  f 
is  performed.  The  frequency  region  of  interest  in  refraction  work  extends  from 
near  zero  to  about  2U  Hz.  Absorption  of  higher  frequency  energy  at  longer 
ranges  sets  this  upper  limit.  The  bandwidth  W  must  be  kept  narrow  so  that 
frequency  selective  phenomena  within  this  20  Hz  band  can  be  discerned.  The 
number  of  significant  Fourier  components  in  a  band,  W,  is  M  =  2WT.  With  the 
RUbt  data,  W  was  3  or  4  Hz  to  maintain  a  reasonable  resolution  in  frequency, 
center  frequencies  were  typically  5,  8,  11,  and  14  Hz.  Since  the  selected  data 

is  modelled  as  a  unidirectional  process  in  white  noise,  in  keeping  with  the 

2irv; .  . 

covariance  expression  in  eq.  4-26,  a  phase  shift: 

is  applied  to  each  matrix  term  to  compensate  for  the  moveout.  In  matrix 


notation,  the  estimated  covariance  matrix  is  expressed  in  two  forms: 

A  I—  W  M  -1 


K 

--  ‘/m  #,  i*  & 

where  H  denotes  the  conjugate  transpose,  w  expresses  the  fact  that  FFT  was 


(4-33b) 


done  on  windowed  data,  K  =  number  of  sensors,  M  =  number  of  frequency 


I 


components  averaged,  and: 


(4-34) 


for  sensor  j,  and: 


... 

for  the  kth  frequency  coefficient  in  W. 

A 

The  term  in  (4-43a)  expresses  the  KxK  S(f)  matrix  as  a  product  of  a 
KxM  and  its  transpose.  The  maximal  rank  of  S(f)  must  then  be  the  lesser  of  M 
or  K.  Thus  if  M*  K,  the  matrix  will  not  be  of  full  rank  and  will  not  be 
invertible  so  that  eq.  4-19  cannot  be  implemented.  Since  narrow  bandwidths 
were  desired,  this  was  the  case  in  all  the  ROSE  experiments  processed.  We  now 
discuss  the  steps  that  were  necessary  so  that  the  MLM  expression  could  still 
be  evaluated  with  S(f)  less  than  full  rank. 

A  third,  algebraic,  expression  for  the  estimated  covariance  matrix 
term  (i,j)  is: 

p.  .  1/2"  \/  /.w'VA  ( 4- 33c ) 

su,  j)  =  '/„■£'  X.  (M)  X  •  Oie.9 

wnere  M  is  the  number  of  frequency  components  and  (i,j)  is  the  phase  shift 
due  to  moveout  at  sensor  i,  frequency  k.  In  the  implementation,  the  matrix  is 


normalized: 


S 


/  .  A  _  _ _ ^  _ _ 

.  .  .  _ _  nufi  i  /  -  c  * _ _ _  1 


(4-36) 


and  the  geometric  mean,  CVAV,  of  the  original  diagonal  elements  calculated  so 


that  the  levels  can  be  restored  afterward: 

C\)A\J -  a) £(i, !)£(;>, a)...  SCk,<) 


(4-37) 


l_c.  H* 
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txperimental ly,  this  scaling  and  normalization  produced  better  results,  making 

A 

Plf,  ✓)  estimates  less  sensitive  to  varying  gains  in  the  sensor  electronics. 

A 

The  diagonal  elements  of  the  normalized  matrix  SWo4tM(i,j),  which  is  still 
singular  if  M<K,  are  unity.  A  pseudo-inverse  of  the  matrix  is  formed  by 
adding  a  small  term,  V  ,  to  all  "ones"  on  the  diagonal  and  then  inverting.  The 
addition  of  this  term  to  the  normalized  matrix  made  it  possible  to  add  the 
same  relative  amount  of  artificial  noise  to  matrices  estimated  from  data  with 
varying  levels  of  energy.  In  practice,  X  ranged  from  .01  to  .04. 

Following  the  calculation  of  steering  vectors,  E ( f ,  <f_x)  for  each 
desired,  the  MLM  expression,  eq,  4-18,  is  evaluated.  After  restoring  levels 
with  CVaV,  the  resulting  estimate  can  be  written  as: 


(4-38) 


p/4  \  __  _ _ CVMV _ 

I 

where  Li(i,j)  is  the  (i,j)  term  of  the  inverse  of  the  matrix  with  (i,j)  term: 


C*#  j.) 


(4-39) 


MLM  Bias 


The  MLM  expression,  eq.  4-19,  produces  an  estimate  of  the  coherent 
energy  across  the  array  (P)  together  with  a  reduction  by  K  of  the  measured 
incoherent  energy  for  the  random  process  upon  which  the  ROSE  data  is  modelled. 
The  estimator  is  well  behaved.  We  have  just  shown,  however,  that  the  actual 
implementation  of  the  MLM  used  does  not  follow  eq.  4-19  exactly.  Since 
estimates  of  energy  partitioning  use  both  the  directional  and  amplitude 
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information  in  P(f,  «/),  it  is  important  that  the  estimated  magnitude  be 
accurate.  A  bias  problem  can  appear  because  of  the  following: 

1)  As  we  have  described,  the  covariance  of  the  process  is  estimated 
from  the  data  itself,  collected  in  a  finite  duration  of  time  and  over  a 
spatially  limited  extent.  We  never  know  the  actual  covariance  of  a  process 
appearing  at  an  array.  The  MLM  estimate  using  this  "covariance"  is  biased. 

2)  The  covariance  matrix  in  many  situations  may  turn  out  to  be 
singular  and  not  invertible,  as  with  the  ROSE  data.  With  the  addition  of 
“white  noise"  terms  on  the  diagonal  of  the  matrix  so  that  a  pseudoinverse  can 
De  formed,  the  behavior  of  the  actual  estimator  does  not  lend  itself  to 
analysis  as  easily  as  eg.  4-19. 


Capon  and  boodman  Bias  Expression 

We  return  to  the  second  expression  for  the  initial  covariance 
estimate,  eq.  4-33b,  before  artificial  diagonal  terms  are  added.  From  a 
probabalistic  viewpoint,  this  equation  is  recognized  as  the  sample  mean 
estimate  of  the  expectation:  , i .e.  the  covariance 

matrix  of  a  K-variate  random  vector  ^  .  Let  the  vector  5,  be  zero  mean, 

complex  uaussian  and  let  sample  vectors  be  normally  distributed  also, 

with  independent  components  for  different  frequencies.  In  this  case,  there  are 
M  independent,  identically  distributed  vector  terms  in  the  sample  mean, 
uoodman  (1963)  shows  that,  under  these  conditions,  the  joint  distribution  of 


all  real  and  imaginary  components  of  MS(f)  has  a  complex  Wishart  density 
function  designated  as:  CW(M,K,  S(f ) ) ,  if  the  matrix  is  of  full  rank  (M>K). 
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For  a  scalar  b  with  density  CW(M,1,1),  the  distribution  is  identical  to  that 
of  a  chi-square  variable  with  the  degree  of  freedom  parameter  equal  to  2M. 
Thus: 

AA  [A-H  for  b  CW(M, 1,1)  (4-40) 


Capon  (ly/U)  shows  that  the  quantities: 


^•=  trX, 
and 


H  P,  CJ.  £) 


is  CW(M, 1,1) 


(4-41) 


£t({x)£ 

HPi  ({.A  . _  is  CW(M-K+1,1,1)  (4-42) 


where  L(f,  )  are  the  steering  vectors  defined  in  Appendix  4-3,  and 
p,  lf.  C  and  »  are  the 


conventional  and  MLM  estimators,  respectively.  Rearranging  and  taking 
expectations,  we  obtain: 


■«t  (4-43) 


for  the  conventional  estimate, and 


-f 


(4-44) 


for  the  MLM  estimate.  Even  if  the  covariance  matrix  has  full  rank,  there  is  a 

M-K+  I 

bias  term:  - — — — ■  in  the  MLM  estimate  due  to  the  fact  that  we  use  a 


M 


windowed  estimate  of  the  actual  covariance.  Although  this  bias  term  is 
tractable,  two  facts  make  it  unsuitable  for  the  present  use,  except  as  a  point 


of  reference: 
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i)  For  M <  K,  the  expression  loses  its  validity.  Even  for  K  =  8 
sensors,  as  in  the  MABS  array,  the  bias  expression  above  is  meaningless,  since 
ivi  :  J  or  4.  Also,  as  we  have  mentioned,  the  resultant  covariance  estimate  is 
not  airectly  invertible. 

ii)  The  bias  term  is  based  on  the  assumption  that  the  terms  in  the 
sample  vectors,  ,  are  normally  distributed  and  independent  in  frequency. 
Appendix  4-b  outlines  a  calculation  of  the  correlation  function  of  two 
components  of 


witn  the  result  that: 


these  vectors: 


1) 


|f,  -fx|>  BW 

jf,  -fj  <■  BW 

f,  -- 

(4-45). 


The  components  X.  (f  )  and  X  (f^)  are  uncorrelated  only  if  they  are  separated 
by  an  interval  larger  than  BW,  the  effective  bandwidth  of  the  window  function 
useq  to  reduce  the  variance  of  our  estimate.  For  a  1  second  cosine  window,  as 
implemented  with  RUSE  data,  BW  is  on  the  order  of  2.5  Hz.,  but  for  T=1  second, 
coefficients  are  spaced  at  about  1  Hz.  Because  of  the  short  data  segments  that 
must  be  used  with  refraction  data,  the  frequency  components  are  therefore 
correlated.  The  bias  expression  due  to  Capon  and  Goodman  is  not  valid.  We  can 
estimate  that  the  effective  number  of  degrees  of  freedom  of  the  chi-square 
variables  b,  and  b^in  Eq.  4-41  and  4-42,  is  reduced  from  the  full  2M  because 
of  tnis,  so  that  the  "M"  term  in  the  bias  expression  is  effectively  decreased. 
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Mooel  Program 

We  nave  shown  that  Capon's  expression  for  MLM  bias  cannot  be  used  with 
the  KObfc  data.  Although  we  can  speculate  on  the  effect  of  correlated 
vectors, ,  no  analytic  results  are  available  for  this  or,  particularly, 
for  the  effect  of  the  artificial  diagonal  terms, V  ,  added  to  the  covariance 
matrix.  To  study  the  influence  on  MLM  bias  by  the  added  artificial  noise,  a 
Monte  Carlo  simulation  was  done.  For  any  given  array  geometry  specified  by  the 
user,  an  MLM  estimate  was  computed  that  was  based  on  the  process  model 
consisting  of  a  sinusoidal  plane  wave  of  amplitude  Jr  and  random  phase, 
together  with  additive,  normally  distributed,  spatially  independent  noise.  The 


Fourier  component  at  sensor  i  and  frequency  k  is  modelled  as: 


(4-46) 


where:  ^  (i)  is  the  random  phase  term  from  a  number  generator  with  a  flat 

distribution  from  0  to  2jr- 

w(i,k)  is  the  complex  “sensor  noise"  term  generated  with  a  zero  mean 
baussian  distribution  with  variance 
<£(i,k)  is  the  moveout  phase  shift. 

The  bandwidth  W  and  number  of  coefficients  (M)  in  the  band  are  chosen  so  that, 
upon  averaging,  the  i,j  term  in  the  simulated  covariance  matrix  is: 

(. »  '/m  X  X*  (4-47) 

0  JBt.s  /  * 

This  expression  is  similar  to  the  corresponding  term  in  the  actual 
implementation  (eq.  4-33c).  In  practice,  M  ranged  from  1  to  76  coefficients 
with  special  attention  given  to  trials  with  M=3  or  4,  numbers  actually  used 


t 


with  the  RUSE  data.  A  constant, o<  is  added  to  all  diagonal  terms  to  replicate 
the  procedure  used  with  real  data.  A  conventional  and  MLM  estimates  for  a 
range  of  is  then  performed  with  the  final  result  expressed  in  decibels: 

lu  log  P(f,  Z^).  The  output  of  each  simulation  is  a  random  variable.  Several 
trials,  usually  ten  or  more,  are  averaged  for  statistical  stability. 

Figures  4-7  through  4-9  show  results  for  a  case  in  which  the  model 
plane  wave  arrives  at  an  angle  of  75  to  a  24  channel,  2400  meter  line  array 
(modelled  after  the  ESP  array).  This  angle  corresponds  to  a  horizontal  phase 
velocity  of  155U  m/sec,  so  that  results  are  for  a  near  endfire  geometry.  The 
nuniDer  of  frequency  components  employed  was  decreased  progressively  from  29  to 
11  to  b  in  a  one  Hz.  band  centered  at  8  Hz.  The  amplitude  P  and  noise 
levels  <5*  and  eC' were  adjusted  so  that  an  unbiased  estimate  at  the  target 
direction  is  15  db,  based  on  the  process  model  and  eq.  4-29.  For  both  the 
conventional  and  MLM  curves,  the  envelope  of  the  mean  estimated  one  sample 
standard  deviation  is  plotted.  In  all  of  these  figures,  at  the  plane  wave 
angle,  the  conventional  estimate  is  considerably  less  biased  than  the  MLM 
estimate  which  decreases  steadily  as  the  number  of  components  used  becomes 
smaller.  However,  at  directions  away  from  target,  the  MLM  estimate  is  sharper 
and  performs  on  the  order  of  5  to  10  db  better  in  sidelobe  rejection 

Since  we  simulated  the  space/time  process  of  a  unidirectional  plane 
wave  in  uncorrelated  (sensor)  noise,  the  theoretical  value  of  the  estimate 
at  'Ct?  *  is  (from  eq  4-31): 

*£*  (4-48) 
where  K  is  the  number  of  sensors.  Diagonal  terms  of  an  estimated  covariance 


t 


array  at  15  off  endfire  (75  ).  Upper  band  is  conventional  and  lower 
band  is  MIM  estimate  of  P.  Unbiased  estimates  would  be  15  db.  Number  of 
frequency  components  averaged  in  this  example  is  29. 


o 
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matrix  are  the  only  locations  in  the  matrix  where  the  variance  of  the  zero 
mean  noise  component  appears.  The  deterministic  oO"  term,  added  to  the 
aiayonal  alone,  is  therefore  included  with  the  actual  variance,  6^  ,  in  eq. 
4-48.  Tnis  expression  proved  to  be  experimentally  correct  since  conventional, 
unDiased  estimates  generated  with  the  simulation  routine  were  generally 
centered  at  P  (f,  ^).The  actual  mean  MLM  estimates  were  always  less  than  this 
value.  Reasoning  that  bias  effects  were  associated  with  random  variables  in 


the  algorithm,  we  model  the  actual  estimate  as: 

P*  +  6l/^l  +  *l/K  (4'49> 

where  represents  an  unknown  bias  coefficient  that  does  not 

effect  the  deterministic  white  noise  term.  We  then  define: 


(4- 50a) 


-to 


Po. 


o>  +  i  (4_50b) 

K  r  U.V  l-*-SN)t  _ 1 

where  5 NR  =  .  The  quantity  A  represents  the  deviation  in  db  of  the 


actual  results  from  the  theoretical,  unbiased  estimate.  Each  time  the  program 

A 

was  run,  a  value  for  A  was  formed.  We  estimate  a  A  from  each  A  by  using: 

A  A  /.«.  .  _ _ _  ' 


{>*«>•***  m.Z'°-u-A 

which  is  a  rearranged  version  of  4-50b.  Figures  4-10  through  4-12  show  a 


(4-51) 


summary  of  the  results  of  this  procedure  for  a  model  line  array  2400  meters 
long  with  8,  12,  and  24  channels  respectively  and  for  at)out  30  off 

oroadside,  corresponding  to  a  horizontal  phase  velocity  of  3000  m/sec.  The 

a  yi/  , 

computed  quantity:  10  log  A  is  plotted  versus  f«t*  .  Each  plot  contains 
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results  for  a  number  of  frequency  coefficients,  M.  The  bias  term  in  Capon's 
formula  loses  validity,  and  the  covariance  matrix  becomes  singular,  when 
M-K+l <  U,  or  when  M<23  for  24,  M<11  for  12,  and  M<7  for  8,  sensors.  We  find 


in  the  empirical  results  that,  for  M  above  these  cutoffs,  the  bias  estimate 

H-K*  I  »  Q 

sot?-  and  that  &  J  kT  =  ipM  at 


approaches  a  constant  value  for  large 


these  values.  Capon's  formula  is  approximately  the  asymptotic  result  as 

oecreases.  below  the  cutoffs,  however,  the  bias  increases  steadily  for 
6-V 

large  .  The  model  program  was  run  for  a  variety  of  different  ratios  of 

SNKs.  For  the  twelve  channel  case  in  fig.  4-10,  the  closely  spaced  groups  of 
points  shown  correspond  to  estimates  at  one  value  of  ,  but  with 

different  SNKs.  In  all  cases,  the  calculated  bias  was  found  to  be  largely  a 
function  of  *  fa  ,  with  relatively  small  sensitivity  to  SNR.  The  curves 
drawn  through  the  points  in  figures  4-10  to  4-12  were  calculated  from  the 
following  expressions: 


uiiuwiuy  eAureibiuns. 

„  *  r-  t  filial 

IDJjv  (3  *  L  i  +cioc*  J 


for  M  >(K-1) 

for  M  £ (K-l) 


(4-52) 

(4-53) 


67,  _ 

where  x  =  /&<  and  £  =  .  These  expressions  were  determined 

after  noticing  the  resemblance  of  figs.  4-10  to  4-12  to  frequency  response 
curves  of  linear  filters.  They  fit  the  bias  data  quite  well  and  were  valuable 
for  simplifying  the  bias  correction  process. 


bias  Corrections 

We  now  describe  the  procedure  followed  in  determining  bias  corrections 


necessary  for  RUSt  data  estimates. 

1)  In  order  to  fit  the  actual  data  to  the  basic  model  used  in  the 
simulation  routine,  we  are  assuming  that,  in  correcting  an  estimated  level  for 
a  significant  event,  the  background  noise  behaves  as  temporally  white, 
spatially  uncorrelated  sensor  noise.  In  order  to  apply  the  empirical  results, 
an  unknown  of  importance  when  working  with  real  data  is  the  value  of  6^  , 
the  power  in  this  uncorrelated  noise  component  in  the  data.  We  evaluate  the 


output  of  the  MLM  program  in  directions  where  no  obvious  coherent  signal  is 

A 

present  and  find  an  average  background  level:  10  log  .  In  this  direction, 
the  original  S(i,j)  matrix,  as  implemented  in  the  MLM  program,  has  estimates 
of  6  on  the  diagonal,  and,  after  normalization,  the!  term  added  to  the 


unity  diagonal  of  (i,j)  is  a  percentage  of  6  .so  that: 


oO-V  6Z 


(4-54) 


Again,  this  expression  for  oC  is  valid  only  in  directions  away  from  any 


significant  propagating  noise. 

y  A 

2)  Letting  'A  =  0  in  this  case,  ^  for  this  ratio  can  be  evaluated 

either  from  figs . 4- 1U  to  4-12  or  from  the  expressions  in  eqs  4-52,  4-53. 

* 

This  ^  is  a  function  of  the  number  of  sensors  and  frequency  coefficients  used 
in  the  estimation  at  the  particular  event  being  corrected. 


3)  For  this  "ambient"  noise  case,  P  =  0,  and  eq.  4-49  becomes: 


(4-55) 


using  eq  4-bb,  we  determine 


"2.  __  (<  Pam& 


(4-56) 
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4)  After  finding  6  ,  we  are  in  a  position  to  calculate  biases  in 

cases  where  P  0.  in  the  actual  implementation,  the  diagonal  terms  of 

i .  j )  become  estimates  of  P  +6  in  the  model,  so  that  the  o*.  terms  in 
tne  model  are  related  to  V  by: 

P*-6X>)  (4-57) 

A  ^  Y 

5)  braphs  of  P  vs.  P  were  desired  for  particular  values  of  s  ,  <*' ,  K, 
and  M  that  were  relevant  to  the  ROSE  data.  For  each  desired  P,  with  Tf  and  6*~ 

1  61/  i  ^ 

being  known,  the  values  of  and  f  a<  is  determined  from  4-57.  A  new  ^  is 

then  evaluated  from  the  graphs  or  the  formulas. 

b)  Using  eq  4-50,  the  biased  P  for  each  P  is  found: 

p  =  S JjP  +  &z4z\ t-  *Yk  (4-58) 

examples  of  final  graphs  of  10  log  P  vs.  10  log  P  are  shown  in  figures  4-13 
and  4-14.  Note  that  for  large  values  of  P,  estimates  were  fairly  unbiased. 

This  is  the  case  since,  for  large  diagonal  terms,  the  effective  constant,  , 

added  to  the  diagonals  is  relatively  large.  The  ratio  ^  /cO'  ,  then,  is 

A 

smaller  and  the  bias  coefficient,  £3  *  is  near  unity.  For  lower  levels,  the 

bias  can  become  quite  large.  A  rough  average  bias  in  the  ROSE  data  for  major 
events  with  energies  significantly  higher  than  the  background  level,  was  found 
to  be  on  the  order  of  5  to  10  db. 

We  have  described  the  fundamental  ideas  behind  the  MIM  algorithm.  To 
illustrate  its  behavior  with  refraction  data,  results  obtained  in  the 
estimation  of  arrival  phase  velocities  from  ROSE  d'.ta  are  presented  in  the 
next  chapter.  We  then  discuss  estimation  of  energy  partitioning  in  Chapter  VI. 


Fig. 4- 13 

Example  of  P  versus  P  curve  (actual  versus  estimated  value  of  coherent 
energy  for  different  values  of  K,  **■  ,  M,  and  Y  . 
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CHAPTER  V 

VELOCITY  ESTIMATION  FROM  ROSE  DATA 

before  discussing  crustal  transmission  loss  estimates  calculated  from 
KUSt  data,  we  first  examine  the  use  of  the  MLM  algorithm  as  an  estimator  of 
crustal  velocities.  This  will  serve  as  an  illustration  of  the  behavior  of  the 
estimator  with  actual  refraction  data.  The  algorithm  resolves  received  energy 
with  respect  to  temporal  frequency  bands  and  horizontal  phase  velocities  (or 
angles  of  arrival)  at  an  array.  If  the  relationship  between  the  horizontal 
phase  velocity  of  coherent  energy  at  the  sensors  and  the  velocity  structure  of 
tne  crust  is  known,  than  the  algorithm  effectively  produces  crustal  velocity 
estimates  as  well.  For  a  horizontally  layered  crustal  model,  the  relationship 
is  quite  simple.  The  layer  sound  speed  will  be  numerically  equal  to  the 
Horizontal  phase  velocity  at  the  array,  which  is: 

~  ~  (5_1) 
where  &  is  the  vertical  angle  of  arrival,  c0  is  the  water  sound  speed,  and  p 
is  the  ray  parameter  or  "slowness"  of  the  arrival.  In  chapter  II  we  showed 
tnat  a  more  realistic  view  of  the  crust  is  based  on  a  model  with  continuous 
velocity  gradients  with  respect  to  depth.  The  horizontal  phase  velocity  of 
coherent  events  in  this  case  is  equal  to  the  sound  speed  in  the  medium  at 
whicn  a  ray  turns  upward. 

Lines  25  and  2L  of  the  ROSE  experiment  took  place  over  thinly 
sedimented  areas  with  a  crustal  age  of  approximately  5  M.Y  (Ewing  & 

Purdy, 1982).  Based  on  a  compilation  of  data  from  529  ocean  basin  refraction 
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experiments,  Christensen  and  Salisbury  (1975)  show  that,  in  relatively  young 
regions  such  as  this,  anomalously  low  mantle  refraction  velocities  (7.1  to  7.8 
km/sec)  are  frequently  observed  and  also  that,  at  offsets  lers  than  35  km, 
"layer  3"  velocities  in  the  range  of  6.7  to  6.9  km/sec  predominate.  First 
arrivals  with  estimated  velocities  in  these  bands  are  the  most  prevalent  in 
the  MLM  analyses.  A  consistent  set  of  second  arrivals,  with  lower  phase 
velocity  estimates,  which  may  be  converted  shear  or,  more  probably,  "layer  2" 
events  are  also  seen  in  line  2S. 

In  analysing  RUSE  velocity  estimates,  two  considerations  must  be  kept 
in  mind,  bathymetry  is  very  complex  near  spreading  centers  and  velocities 
estimated  from  array  data  are  influenced  by  topography  to  the  extent  that, 
without  appropriate  corrections,  errors  on  the  order  of  1  km/sec  may  occur, 
bathymetric  data  with  sufficient  resolution  for  correcting  this  problem  is  not 
available  for  the  ROSE  experiment  because  of  navigation  failure.  Secondly,  ESP 
lines  2b  and  2L  crossed  at  least  two  fracture  zones  (see  fig.  5-1  from  Purdy 
(1982)).  Results  obtained  by  Purdy  (1982)  from  OBS  data  near  these  fracture 
zones  are  compatible  with  an  increased  thickness  of  low  velocity  material  in 
the  uppermost  crust,  a  feature  of  fracture  zone  troughs  (Ludwig  &  Rabinowitz, 
(lyau);  Detrick  and  Purdy, ( 1980) .  This  non-homogeneity  in  the  structure  must 
be  taken  into  account  in  viewing  the  MLM  results  from  the  standpoint  of 
"normal"  crustal  models. 

In  this  chapter,  following  illustrations  of  some  time  profiles  of 
horizontal  and  vertical  array  data,  contoured  plots  of  the  relative  strength 
of  arrivals  with  respect  to  phase  velocity  and  travel  time  are  presented.  A 
summary  plot  of  all  experiments  in  one  shooting  line  is  then  discussed  with 


special  attention  given  to  arrivals  beyond  the  first.  Despite  the  effects  of 
rough  topography  mentioned  above,  we  show  that  the  ability  of  the  array  to 
discriminate  arrivals  with  different  relative  phase  velocity  can  still  be 
applied.  A  travel  time-offset  plot  with  events  labelled  with  respect  to  their 
approximate  phase  velocity  range  is  presented  for  both  ESP  lines. 

Time  Prof i les 

He  first  look  at  some  of  the  raw  data  as  it  was  entered  into  the 
velocity  analysis  program.  In  figures  5-2a  and  2b,  12  of  the  24  channels  of 
data,  arranged  sequentially,  from  the  horizontal  array  (ESP)  are  plotted 
versus  travel  time  (vertical  axis).  Horizontal  offset  for  this  data  was  26  km. 
and  each  tick  mark  represents  one  second.  The  first  arrival,  at  about  two 
seconds  into  the  record,  appears  almost  simultaneously  on  all  the  sensors, 
i.e.  with  little  "moveout".  Since  it  appears  coherently  across  the  array,  this 
high  phase  velocity  event  is  a  "refracted"  arrival  (IP)  emitted  from  the 
seabed,  meeting  the  array  almost  broadside  (i.e.  at  a  vertical  angle  close  to 
zero,  or  a  "grazing"  angle  near  90*).  The  next  visually  apparent  arrival 
occurs  about  4.6  seconds  later  and  displays  the  same  small  moveout  across  the 
elements.  This  is  a  "refracted/  reflected"  arrival  (2P),  as  defined  in  Chapter 
II.  The  water  depth  in  this  area  is  about  3  km  so  that  the  round  trip  time 
from  the  seabed  to  the  water  surface  and  back  is  on  the  order  of  4  seconds. 
Turning  to  figure  5-2b,  two  stronger  arrivals  occur  at  T=17  and  T=19  on 
channel  1,  the  sensor  closest  to  the  receiving  ship  and  furthest  from  the 
shot.  Unlike  the  previous  two,  these  are  displaced  in  time  across  the  array. 


\Z  cnannels  of  data  from  the  ESP  array  plotted  versus  travel  time.  Shaded 
areas  in  this  figure  and  the  next  indicate  typical  data  windows  used  when 
frequency/wavenumber  function  estimates  are  made  for  high  and  low  phase 

velocities  respectively. 
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i.e.  with  larger  moveout.  They  arrive  at  a  smaller  grazing  angle,  i.e.  more 
from  the  shot  direction  than  the  seafloor.  These  are  water  arrivals:  1W  and 
2to,  with  low  horizontal  phase  velocities. 

Figure  5-3  shows  data  from  the  12  channel  vertical  array  (MABS)  at  a 
range  of  17  km.  As  we  have  discussed,  channels  5,  7,  11,  and  12  were 
malfunctioning.  Unlike  the  ESP  case  above,  the  first  arrivals  here,  at  T=6.3 
sec  on  the  shallowest  channel,  #1,  have  the  greatest  moveout.  This  is  expected 
for  energy  coming  from  the  sea  floor  direction  and  arriving  almost  endfire. 
Another  set  of  arrivals  appears  at  T=  7.8  seconds  on  channel  1.  This  event 
also  has  a  large  moveout,  but  it  propagates  in  the  opposite  direction.  The 
shallowest  channel  in  the  deployment  was  at  a  depth  of  about  1  kilometer.  If 
we  were  to  "continue"  the  locus  of  first  points  at  each  sensor  for  both  of 
these  events  up  to  a  hypothetical  sensor  at  the  surface,  as  indicated  by  the 
dashed  lines  in  the  figure,  the  lines  intersect.  Both  arrivals  result  from  the 
same  crustal  "refraction"  (IP)  with  the  later  one  being  caused  by  reflection 
at  the  surface.  At  about  9  seconds  into  the  data  on  channel  1,  another  endfire 
set  of  arrivals  can  be  seen  across  the  sensors.  Occurring  at  a  4.5  second 
interval  after  IP,  this  is  again  a  refracted/ref lected  2P.  Finally,  a  set  of 
arrivals  that  appear  almost  horizontally  across  the  data  is  seen  at  T  =11.2 
seconds.  Because  of  their  large  amplitudes,  the  tape  recorders  saturated  at 
this  point,  and  there  is  no  usable  information  beyond.  Since  there  is  no 
moveout  for  this  set,  however,  it  is  clearly  a  direct  water  arrival. 


KtSULTS  OF  VELOCITY  ANALYSIS  PROGRAM 

toe  now  present  some  results  of  the  MLM  procedure  applied  on  the  data 
just  discussed.  Figure  5-4  is  a  contoured  plot  of  the  estimated  wavenumber 
function,  in  a  three  Hertz  band  centered  at  8  Hz,  versus  phase  velocity 
vertically)  and  travel  time  (horizontal).  These  estimates  were  calculated 
from  the  horizontal  array  data  shown  in  figs.  5-2.  At  each  time  T,  the  target 
angle,  looking  downward,  was  stepped  over  a  range  of  0*  to  about  80*  (grazing) 
in  increments  which  correspond  to  equal  slowness  (p)  intervals  of  about  5.8 
yts/meter.  The  range  of  phase  velocities  is  1.5  to  8.8  km/sec.  The  contouring 
of  the  levels  was  done  at  2  db.  intervals.  Proceeding  from  left  to  right,  a 
background  level  of  -50  db.  quickly  changes  to  a  sharp  peak  about  two  seconds 
into  the  figure.  This  peak,  at  a  phase  velocity  of  6.8  km/sec  corresponds  to 
the  first  arrival  (IP)  observed  on  the  profile  in  fig.  5-2.  The  maximum  level 
here  is  at  about  -15.5  db.  The  value  of  the  level  estimates  has  not  been 
corrected  for  the  effects  of  MLM  bias.  A  second  event  ("  IP*  "),  not  visually 
apparent  in  fig.  5-2,  occurs  a  fraction  of  a  second  later  at  a  slightly  lower 
Cp.  Tnis  "doublet"  phenomenon  is  seen  frequently  in  the  ROSE  data.  One 
possible  explanation  is  the  existence  of  a  low  velocity  zone  at  the  base  of 

the  crust  (Lewis  and  Snydsman,  1977).  Evidence  of  a  low  velocity  region  in 

lower  crust,  from  an  OBS  experiment  conducted  near  the  East  Pacific  Rise,  is 
also  presented  by  Orcutt  (1976). 

The  next  prominent  arrival,  at  T  =11.2  seconds,  was  also  not 
discernable  on  the  time  profile.  With  a  phase  velocity  of  about  4.5  km/sec, 

and  a  level  about  13  db  lower  than  IP  and  IP',  it  is  either  a  late  arrival 
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Fiq.5-4a 

Results  of  velocity  analysis  of  data  in  figs.  5-2.  Horizontal  axis 
is  travel  time.  Vertical  axes  are:  ljapparent  phase  velocity  across 
array,  2)equivalent  slowness  (p)  of  arrival,  and  3)  equivalent  grazing 
angle  at  array  for  plane  waves.  Amplitude  of  estimates  are  contoured 

in  4  db  increments. 
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See  caption  for Fig.5-4a. 


from  an  area  with  a  lower  sound  speed  ("layer  2",  if  we  use  the  Raitt  model), 
or  a  converted  shear  wave  (IS),  although  the  phase  speed  is  somewhat  high  for 
the  latter  case.  The  next  peak,  4.5  seconds  after  IP,  is  the 
refracted/ ref lected  2P  with  a  level  of  about  -18  db  and  a  velocity  of  7.3 
km/sec.  The  level  of  this  arrival  here  is  not  higher  than  that  of  IP,  although 
this  is  often  the  case.  An  echo  of  IP'  appears  next  at  T= 1 2 . 7  followed  a 
second  later  by  a  weaker  (-27  db)  event  at  6.2  km/sec.  Again,  these  last  two, 
although  considerably  stronger  than  the  background  level,  were  nevertheless 
not  visually  discernable.  Finally,  a  progression  of  arrivals  begins  after  T= 1 7 
seconds  at  very  shallow  grazing  angles.  The  fact  that  each  of  these 
progressively  increases  in  angle  is  in  accord  with  the  interpretation  of  these 
as  water  bounces  with  higher  order  reflections  encountering  the  array  at 
larger  angles. 

Before  turning  to  MLM  estimates  of  the  MABS  data,  we  first  present 
results  from  a  Line  2S  shot  at  a  31  km  offset  in  fig.  5-5.  The  resolution  in 
phase  velocity  of  prominent  events  is  not  as  great  as  in  the  preceeding  case. 
This  is  due  to  the  fact  that  only  12  of  the  24  data  channels  were  used  when 
processing  line  2S.  For  this  experiment,  the  doublet  phenomenon  is  again  seen 
at  T=b.2  and  T=9  seconds,  as  well  as  the  2P  arrival  4.5  seconds  later.  With 
the  increased  offset,  a  "3P"  arrival  also  occurs  at  T=17.5  followed  by  3  water 
arrivals  at  T=20,  21,  and  24.5. 

In  Figs.  5-6a  and  5-6b,  we  look  at  MLM  estimates  obtained  with  the 
vertical  (MABS)  data  of  fig.  5-3.  This  data  was  obtained  from  the  same  shot  as 
the  ESP  data  just  discussed,  although  the  offset  to  the  MABS  was  only  17  km. 
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Again,  vertical  resolution  (in  phase  velocity)  is  not  as  sharp  as  in  the  ESP 
2l  results  due  to  the  smaller  number  (8)  of  sensors,  and  the  fact  that,  for 
high  phase  velocity  arrivals,  the  energy  encounters  the  array  close  to 
endfire  The  IP  arrival  at  6.2  seconds  on  channel  1  of  fig.  5-3  is  represented 
as  the  peak  at  T  =6.8  sec.  on  the  contour  plot.  Since  the  mathematical  origin 
of  the  array  geometry  used  in  the  algorithm  is  at  the  surface,  for  a  channel 
at  depth  z^  and  a  target  angle  looking  below  the  vertical  array,  an  estimate 
at  time  t  uses  data  that  appeared  on  channel  i  at: 

jt  ~  (5-2) 

'-o 

For  channel  1,  at  1  km  depth,  this  is  .66  sec  at  a  phase  velocity  of  7000 
km/sec.  At  T= 1 1  seconds,  the  event  at  phase  velocity  5  km/sec,  with  amplitude 
20.6  db,  is  the  upgoing  2P  arrival  observed  in  fig  5-3.  In  this  case,  the  2P 
level  is  indeed  higher  than  that  of  the  IP.  In  fig  5-3,  we  saw  that  beginning 
at  T»ll  seconds,  the  recorded  MABS  data  was  not  usable  for  velocity  analysis. 
The  2P  event  estimated  at  this  time,  however,  is  still  based  on  earlier, 
coherent  data.  Using  eq  5-2,  the  approximate  locus  of  invalid  data  is 
indicated  in  the  shaded  region. 

Although  the  horizontal  offset  in  this  data  is,  respectively,  9  and  14 
km  less  than  the  offset  of  the  2L  and  2S  data  discussed  above,  the  higher 
level  (+17.3  db)  of  IP  in  these  estimates  is  partly  due  to  a  decreased  MLM 
bias  for  the  smaller  number  of  channels  used.  In  the  next  chapter, 
compensations  for  this  effect  are  calculated  so  that  levels  are  made 
insensitive  to  channel  number  for  the  8,  12,  and  24  sensor  data. 

In  figure  5-6b  the  target  angle  used  in  the  MLM  alcorithm  was  set  for 
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Results  of  velocity  analysis  on 8  channel  vertic< 
at  angle  looking  down.  Shaded  area  is  approximate 
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estimation  in  the  same  phase  velocity  range  as  above,  except  that  the  array 
looks  upward.  Arrivals  from  the  surface  are  contoured.  The  IP  reflected 
arrival  appears  earlier  in  this  case  since  data  subsequent  to  the  estimation 
time  is  used.  The  locus  of  invalid  data  has  increased  in  this  case,  as 
indicated  in  the  shaded  region  The  peak  level  of  IP  here  is  about  3db  higher 
than  that  of  the  corresponding  event  in  5-6a.,  due  to  interference  effects  at 
tne  surface  and  statistical  fluctuations. 

For  an  overall  view  of  the  behavior  of  the  velocity  estimates 
generated  by  the  routine  at  different  offsets,  figure  5-7a  and  5- 7b  are 
schematic  outlines  of  the  main  events  from  contour  plots  of  ten  shots  in  one 
ESP  refraction  line  (2S)  at  the  8  Hertz  center  frequency.  Each  of  the  ten 
bands  gives  a  summary  of  one  experiment  with  the  prominent  events  plotted 
horizontally  with  travel  time  after  the  first  arrival,  and  vertically  with 
estimated  phase  velocity  from  zero  to  ten  km/sec.  The  estimated  levels  and 
path  designations  (where  possible)  are  annotated  at  each  peak.  All  ten  shots 
appear  with  first  arrivals  alligning  vertically.  The  appearance  of  lltf  in  each 
band  (indicated  by  the  dashed  lines)  occurs  at  the  first  event  with  relatively 
low  ( l.b  to  2  km/sec)  estimated  velocity.  Preceeding  1W,  most  of  the  energy  is 
seen  to  concentrate  at  time  intervals  of  4.5  seconds  and  nh  estimated 
velocities,  due  to  refraction/reflection.  There  is  a  "spli  of  IP  into 

distinct  multipaths  in  most  experiments  and  this  phenomemon  is  also  seen  in 
tne  2P  regions.  "Medium"  velocity  arrivals  (in  a  4  to  6  km/sec  band)  occur  two 
to  three  seconds  after  IP  in  four  of  the  bands.  The  level  of  peaks  in  the  2P 
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regions  is  not  consistently  higher  or  lower  than  the  IP  levels.  This  could  be 
due  to  interference,  varying  from  constructive  to  destructive,  due  to 
bathymetric  variations  which  were  on  the  order  of  a  wavelength  (200  meters). 
After  the  water  wave  arrivals,  most  shots  exhibit  the  multiple  structure  we 
have  discussed,  gradually  increasing  in  "velocity"  and  separation  from  the 
preceding  multiple. 

Within  each  shot  band,  the  estimates  give  a  good  indication  of  the 
relative  variation  in  phase  velocity  of  coherent  arrivals.  This  has  made  it 
possible  to  identify  different  types  of  arrivals  in  one  experiment.  Actual 
numerical  velocity  estimates  from  an  array,  however,  are  suspect  in  regions 
where  the  horizontal  layer  model  is  not  valid  due  to  rough  topography.  For 
instance,  if  we  take  the  simple  one  inteface  model  and  allow  the  boundary  to 
be  slightly  inclined  from  the  horizontal  by  ,  the  variation  of  estimated 
phase  velocity  in  the  water  with  the  inclination,  via  eq.  5-1,  is: 

ACp  -  Cp  (5-3) 

The  sensitivity  for  this  simple  model  is  large  for  events  with  high  phase 
velocities.  At  c^  =”  7  km/sec,  for  instance,  an  inclination  of  one  degree  would 
change  the  estimated  velocity  by  500  m/sec.  If,  in  an  expanded  model,  a  gently 
varying  layer  at  the  seabed  was  above  a  set  of  strictly  horizontal  layers, 
than  phase  velocity  estimates  at  the  array  could  be  corrected  by  adequately 
sampling  the  bathymetry  along  the  line  and  using  eq.  5-3  to  correct  for 
various  inclinations  encountered.  Because  of  the  proximity  of  the  experiment 
to  the  East  Pacific  Rise,  however,  the  bathymetry  in  the  ROSE  area  was 
extremely  complex.  Figure  5-8  is  a  diagram  of  sampled  depths  for  line  2S,  at 
the  locations  of  the  emerging  rays,  with  approximate  bottom  inclinations. 
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oottoni  slopes  in  the  2  to  3  range  were  not  uncommon  and  phase  velocity 
estimates  with  the  simple  model  just  mentioned,  could  be  more  than  1  km/sec  in 
error.  Unfortunately,  we  obtained  depth  samples  only  at  intervals  of  1  km  or 
more  (on  the  order  of  the  ESP  array  length),  and  attempts  to  correct  the 
velocity  estimates  have  not  been  fruitful. 

In  this  situation,  although  raw  velocity  estimates  from  the  velocity 
spectral  analysis  routine  are  suspect,  the  ability  to  discern  the  relative 
difference  in  phase  velocities  of  sequential  events  is  still  of  use.  In 
figures  5-9a  and  5-9b,  travel  time/offset  plots  for  two  ESP  lines  (2S  and  2L) 
are  shown.  These  figures  were  constructed  from  range  information  generated  by 
the  RaYuIST  unit  and  travel  times  based  on  the  first  water  arrival  ( 1W ) . 
Altnougn  a  more  accurate  system  for  measuring  arrival  times  would  be  required 
for  in  depth  analysis,  the  ability  to  discriminate  relative  phase  velocity  at 
tne  array  did  produce  fruitful  results.  Estimated  velocities  are  divided  into 
three  categories:  low  (1500  to  3500  km/sec),  medium  (3500  to  5500),  and  high 
(5500  and  above).  The  suite  of  prominent  arrivals  in  time  are  plotted 
vertically.  The  "doublet"  (and  sometimes  "triplet")  phenomenon  of  closely 
spaced  high  velocity  events  are  indicated  in  the  circles.  In  fig.  5-9a,  with 
fairly  dense  shot  spacing,  we  were  able  to  discern  two  distinct  first  arrival 
slopes.  The  first  at  6.6  km/sec  would  correspond  to  the  approximate  sound 
speed  in  layer  3  while  the  slope  (7.8  km/sec)  at  the  largest  offset  indicates 
a  MOHu  refraction.  Since  upper  mantle  velocity  is  normally  8.2  km/sec  (Lewis, 
1978)  in  older  regions,  this  lower  estimate  is  in  accord  with  the  "anomalous 
mantle"  in  fig.  23  of  Christensen  &  Salisbury  (1975). 


fig  5- 9b 


•  HIGH  VELOCITY  IP,  2P  TYPE  ARRIVALS 

(•)  MEDIUM  VELOCITY,  LOWER  AMPLITUDE,  (POSSIBLE  SHEAR)  ARRIVALS 
X  VERY  LOW  PHASE  VELOCITY  (WATER  WAVE)  ARRIVALS 


Fig,5-9b 

Travel  time/  offset  diagram  obtained  for  line  2L. 
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A  consistent  set  of  "medium1’  MLM  estimates,  occurring  after  the  first 
arrival  with  an  approximate  slope  of  5.7  km/sec,  confirms  the  presence  of 
"layer  2"  events  .  Because  average  crustal  shear  velocity  is  about  4  km/sec, 
the  possibility  of  these  "medium"  velocity  events  being  converted  shear 
arrivals  is  ruled  out.  The  refraction/ref lection  (2P)  events  are  seen  in  the 
line  parallel  to  the  locus  of  first  arrivals.  As  many  as  four  low  phase 
velocity  water  arrivals  can  also  be  seen,  the  lines  formed  by  these  events  on 
the  T-X  diagram  all  tending  to  converge  at  large  offsets. 

In  figure  5-9t>,  although  the  refraction  line  was  actually  run  out  to 
ranges  in  excess  of  UJO  km,  water  wave  data  was  not  available  to  us  beyond  80 
km,  at  which  point  travel  time  calculations  could  not  be  continued.  First 
arrivals  indicate  both  a  6.6  km/sec,  and  a  higher  (8.9  km/sec)  slope, 
intersecting  at  a  range  of  30  km.  This  extremely  high  mantle  velocity  estimate 
is  due  to  errors  caused  by  calculating  first  arrival  slopes  from  sparsely 
sampled  data.  In  this  line  there  is  only  one  "medium"  velocity  event,  at  26 
km,  tnat  may  be  a  layer  2  arrival.  The  "doublet"  phenomena  is  especially 
prominent  in  2P  and  3P  refract  ion/ref  lections,  and,  at  40  &  52  km,  a  "4P" 
arrivals  occur. 


Although  extremely  rough  bathymetry  and  errors  in  position  information 
lessen  the  accuracy  of  velocity  estimates  from  one  shot  alone,  we  have  shown 
that  knowledge  of  the  relative  values  of  velocity  estimates  in  one  experiment 
can  still  aid  in  the  interpretation  of  a  travel  time-offset  diagram  of 
refraction  data.  Furthermore,  we  have  been  able  to  identify  events  beyond  the 
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first  arrival  and  can  discriminate  arrivals  from  different  trajectories  that 
appear  simultaneously  at  a  receiver. 

The  interpretations  of  the  MLM  estimates  of  the  first  arrivals  in  the 
T-X  plot  are  supported  by  the  fact  that,  at  a  range  of  30  to  40  km,  the 
replacement  of  layer  2  events  as  first  arrivals  by  higher  velocity  mantle 
refractions  is  a  common  occurrence.  An  increase  in  amplitude  at  this  distance 
occurs  frequently,  as  discussed  in  Chapter  II.  In  chapter  VI,  we  use  level 
estimates  from  the  MLM  routine,  calculate  crustal  transmission  losses,  and 
determine  estimates  of  this  energy  focusing  at  these  ranges. 
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CHAPTER  VI 

TRANSMISSION  LOSS  CALCULATIONS 

A 

We  now  describe  the  procedure  followed  In  talcing  values  of  Pff,/1) 
generated  by  the  MLM  algorithm  and  calculating  numerical  estimates  of 
transmission  loss  for  ROSE  refraction  events.  We  have  shown  that  the 
frequency /wavenumber  estimates  represent  energy  arriving  at  the  array 
partitioned  with  respect  to  both  temporal  and  spatial  frequencies  (power  per 
Hertz  per  steradian).  The  relation  between  the  estimates  and  the  acoustic 
quantities  defining  transmission  loss  is  first  presented.  We  then  discuss  five 

A 

corrections  that  are  applied  to  P(f,£)  so  that  valid  TL  estimates  are 
produced.  The  method  followed  in  the  calculation  of  source  level  (SL)  for  each 
shot  is  also  described.  Estimates  of  transmission  loss  versus  range  are  then 
presented  for  ESP  lines  2S  and  2L.  The  particular  paths  for  which  transmission 
losses  are  calculated  are  IP,  2P,  1W,  and  the  "layer  2"  arrivals. 


Relation  of  Transmission  Loss  and  Frequency /Wavenumber  Estimates 

Transmission  loss  at  a  point  r  is  defined  as  (from  Clay  A  Medw1n,1977) 


TL  C-i')  -  SL  -  SPl  (ji\ 


(6-11 


SL  denotes  the  source  level: 


SL- l0^ 


(6-2) 


where  p%  (R^  )  is  the  mean  square  pressure  at  reference  distance  R,^ 
SPL(r)  Is  the  sound  pressure  level  at  r: 


(6-3) 
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In  tnis  chapter,  the  reference  distance  used  is  R„.,  =  1  meter,  and  the  root 


mean 


square  reference  pressure  is  p  =  1  ^Pa.  The  mean  square  pressure  is 
related  to  the  spectral  density,  ^^(f,r)  by: 

(6.41 

If  the  sound  field  is  modelled  as  a  stationary  random  process,  then  Jfjf,  r) 
is  equivalent  to  the  spectral  covariance  function  of  the  random  pressure 
process:  (f,r)  =  Sp(f,  ci-r^)  at  rA  =  ri  =  r.  Both  describe  the  density 

of  energy  with  respect  to  temporal  frequency,  f,  at  location  r. 

For  a  homogeneous  process,  the  covariance  function  can  be  written  as: 


(6-5) 


where  P^lf,  *)  is  the  frequency/wavenumber  function  of  the  random  pressure 
process.  The  mean  square  pressure  can  then  be  written  as: 


(6-6) 


pp(f»'C)  thus  represents  the  density  of  energy  per  Hz.  per  steradian. 

In  Chapter  IV,  an  acoustic  event  resulting  from  an  explosive  source  in 

a  refraction  experiment  was  modelled  as  a  windowed  segment  of  a  unidirectional 

plane  wave  pc(*),  propagating  in  a  direction  >fp/  1*1-  The  frequency/ 

wavenumber  function  of  this  model  process  is  impulsive,  i.e.  Pp(f,vT_)  = 

P0(f)  £('C">£p)*  where  Pe(f)  is  the  power  spectrum  of  P0  ( * ) -  By  substitution, 

eq.  b-o  can  then  be  rewritten  as: 

<Xi 


w  <30 

'1*1  fa  fa) 


(6-7) 
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ln  the  MLM  algorithm,  estimates  of  the  frequency/wavenumber  function  are  made 
for  a  set  of  center  frequencies,  fc  ,  in  bands  of  width  W.  We  assume  the 
estimated  function  is  approximately  constant  over  W  so  that: 


*  p.  ({ .O  w 


(6-8) 


'p^w(r)  is  the  mean  square  pressure  in  a  frequency  band  centered  at  f0  with 
width  W. 

In  estimating  transmission  loss,  we  choose  an  event  represented  by  a 

A 

large  value  of  the  estimated  wavenumber  function,  P(f,y$)  and  consider  the 
event  as  an  arrival  of  the  model  process  at  \Qt=  The  transmission  loss  for 
the  chosen  event  at  is  then: 


(6-9) 


i  L-^#vo (<£)  ~  £  L  -  S  PL (-&) 


v -  \  f  fV 

(6-10) 

We  choose  *f^  =  W  so  that:  *  A 

This  is  a  modified  version  of  eq  6-1  with  TL^^y  (r)  representing  the  loss  in 
a  specific  frequency  band  arriving  at  the  receivers  at  the  angle 
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Corrections  to  MLM  Estimates 

before  eq  6-11  can  be  implemented,  we  must  relate  P^(f,  tft)  ,  the 
frequency/wavenumber  function  of  the  model  pressure  field  with  P^(f ,  the 
quantity  actually  estimated  in  the  MLM  algorithm  from  the  sampled  data  points 
on  the  magnetic  tapes.  Five  corrections  to  P^f.y)  are  needed  to  remove 
effects  of:  i)  sampling  in  time  and  frequency;  ii)artifacts  of  receiver 
location  (Lloyd  mirror  effect);  iii)  array  gain,  iv)  hydrophone  sensitivity, 
and  v)  MLM  bias. 


tt6n  dS:  / 


Sampling  Correction 

In  Appendix  6-1,  it  is  shown  that  the  spectral  density  function,  for 
windowed  segments  of  length  T,  sampled  at  the  Nyquist  interval  *  T,  can  be 
written  as:  ..  l"***.,  'Hr 

%r  h  J  l6-’2> 

M  is  the’  number  of  Fourier  coefficients,  P(kAf),  in  W.  k„  is  the  coefficient 
number  corresponding  to  fe .  If  we  compare  eq  6-12  with  eq  4- 34c,  the  term  in 
the  above  bracket  is  recognized  as  the  implemented  expression  for  the  diagonal 
elements  of  the  estimated  covariance  matrix,  rvn  .  In  order  that  the 
estimated  matrix  be  equivalent  to  the  spectral  density  function,^  (f,r) ,  a 
correction  due  to  sampling:  e=  AT^/T,  is  applied  to  the  matrix.  This  same 
correction  must  be  applied  to  the  MLM  frequency/wavenumber  estimates. 

In  the  ROSE  data  set,  T  =  .004  sec.,  and  the  effective  data  length 

(after  windowing)  was  T~.b  sec.  so  that  a  correction  of: 

(.oo*+)x 


1 0 


io  JUy 
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was  subtracted  from  10  log  P(f,  )  for  al_^  of  the  experiments. 


Lloyd  Mirror  Correction 

The  Lloyd  mirror  effect  is  important  in  characterizing  the  sound  field 
near  a  free  surface,  especially  at  low  frequencies,  where  wavelengths  are 
greater  than  the  dimension  of  the  average  surface  roughness  so  that  reflection 
is  specular.  A  correction  for  this  effect  was  necessary  for  data  received  with 
tne  horizontal  (ESP)  array.  At  an  approximate  depth  of  10  meters,  surface 
reflections,  when  added  to  arrivals  from  below,  significantly  alter  the 
amplitude  of  the  waveform  at  a  sensor.  Referring  to  Fig.  6-1,  EDA  represents  a 
pure  sinusoidal  plane  wavefront  with  frequency  f,  arriving  at  a  sensor  at 
point  A  at  time  t.  The  vertical  angle  is  and  hydrophone  depth  is  BA  =  d. 

The  sensor  will  also  be  influenced  at  this  moment  by  a  surface  reflected 
arrival  that  has  traversed  the  extra  distance  DC  +  CA.  The  necessary 


correction  for  this  effect,  as  shown  in  Appendix  6-2,  is: 


aoi»  |; 


7ir 


(6-13) 


Typical  correction  curves  for  different  frequencies,  f,  versus  vertical 
angle,  oC  ,are  shown  in  fig.  6-2.  Since  estimates  were  performed  across  a  3  or 
4  hz  oand,  corrections  were  averaged  across  each  band.  For  high  phase 
velocities,  this  factor  is  on  the  order  of  3  db  at  8  Hz.  For  each  event,  the 
estimated  angle  of  arrival  obtained  with  the  MLM  algorithm  was  used  in  eq. 
b-13.  The  frequency  averaging  is  particularly  important  at  large  vertical 
angles  (water  arrivals)  since  corrections  are  large  and  sensitive  to  frequency 


where  predicted  amplitudes  are  small.  Because  of  this  severe  attenuation  near 
endfire,  the  bottom  reflected  water  wave,  rather  than  the  essentially 
unobservable,  higher  angle  direct  arrival,  is  observed  in  ESP  data. 


Array  Directivity 

Another  correction  applicable  to  the  ESP  array  only  is  due  to  the  fact 
that  each  of  the  24  "sensors"  was  actually  composed  of  two  fifty  meter 
streamer  sections  connected  in  parallel.  To  correct  for  directivity  effects  of 
each  of  these  small  "arrays",  each  channel  is  modelled  as  an  unphased  100 
meter  long  line  array  with  beam  pattern,  from  eq.  4-12: 

<6-14) 

where  ^  is  the  vertical  angle.  The  correction  applied  to  the  data  due  to 
this  effect  is:  _ , 

0.0  Ap  L<^  J  (6-15) 

which  is  plotted  in  fig.  6-3.  As  with  the  mirror  effect,  high  velocity 
arrivals  at  low  angles  are  not  heavily  affected  by  this  correction,  but  they 
are  quite  sensitive  for  water  arrivals.  Frequency  averaging  across  the 
estimation  band  was  done  for  all  events. 

hydrophone  Sensitivity 

The  data  recorded  on  the  ESP  and  MABS  tapes  are  the  voltages  that  were 
present  at  the  output  of  the  acquisition  systems.  Hydrophone  sensitivity 
corrections  are  necessary  to  convert  this  data  into  units  of  sound  pressure. 

Figure  6-4  is  a  schematic  illustration  of  the  results  of  a  calibration 


Fig.  6-3 

Array  directivity  effect  for  100  meter  line  array  at  3  frequencies 
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TABLE  6-1 


EXPERIMENT 

ARRAY 

K 

CENTER  FREQ 

M 

* 

2S 

MABS 

8 

8 

4 

.03 

12 

4 

.03 

16 

4 

.03 

2S 

ESP 

12 

8 

4 

.02 

12 

4 

.02 

16 

4 

.02 

2L 

ESP 

24 

5 

3 

.03 

(to  45km) 

8 

3 

.03 

11 

3 

.03 

14 

3 

.03 

2L 

ESP 

24 

5 

3 

.04 

(at  52.5  km) 

8 

3 

.04 

11 

3 

.04 

14 

3 

.04 

2L 

ESP 

24 

5 

3 

.01 

(52.5  to  104km) 

8 

3 

.01 

11 

3 

.01 

14 

3 

.01 
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performed  on  the  MABS  array  prior  to  deployment.  By  averaging  over  the  eight 
functioning  channels,  sensitivities  of  about  -125  and  -121  db  re  1  volt/j^Pa 
were  found  for  the  8  Hz.,  and  the  12  &  16  Hz  bands,  respectively. 

The  sensitivity  of  the  phones  of  the  ESP  array  were  given  as 
-18b  do  re  IV  per^Pa.  For  impedance  matching  purposes,  the  coupling  of  each 
streamer  section  to  the  ship  made  use  of  a  9:1  (18  db)  step  down  transformer, 
so  that  the  effective  sensitivity  was  about  -203  db.  Values  appearing  on  the 
final  tapes  sent  to  Woods  Hole  are  in  millivolts  so  that,  for  data  processing 
purposes,  the  effective  sensitivity  for  the  ESP  array  is: 

1-ZU3  +  ou)  =  -143  db  re  IV/ycPa. 


MlM  bias 


Table  6-1  outlines  the  important  parameters  that  were  used  in  the  MLM 

processing  and  bias  calculation  for  the  ROSE  data.  The  "white  noise" 

factor,  V  was  varied  during  the  processing  of  line  2L. 

Using  the  correction  procedure  discussed  in  Chapter  IV,  we  first  used 

z 

the  background  levels  of  each  shot  to  provide  us  with  an  estimate  of  6  .  This 
was  done  for  all  experiments  and  in  all  the  various  frequency  bands.  The 
results  were  then  averaged  to  provide  the  following  values  of  6 *  used  for 
bias  correction: 


Line  2L,  ESP  array:  6  =  .004  at  5  Hz 

•  =  .0035  at  8  Hz 

"  =  .003  at  11  Hz 

"  «  .002  at  14  Hz 

Line  2S,  ESP  array:  6*  =  .0035  at  8  Hz 

*  »  .003  at  12  Hz 

•'  =  .002  at  16  Hz 


NJO. 
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Line  2S,  MABS  array:  6  =  .175  at  8  H2 

••  =  .175  at  12  Hz 

<«  =  .175  at  16  Hz 

The  difference  in  levels  between  the  MABS  and  ESP  is  about  17  db 


corresponding  to  the  difference  in  hydrophone  sensitivities.  With  6  known,  P 

v\ 

versus  P^^^curves  {as  in  figures  4-15  and  4-16)  were  drawn  for  all  necessary 
combinations  of  V  ,  M,  6*  ,  and  K  needed  in  Table  6-1.  and  the  necessary 


corrections  were  found  from  these  curves. 

Since  MLM  bias  is  a  function  of  ^  ,  a  test  of  the  bias  correction 
procedure  was  performed  with  actual  data  using  different  values  of  this 
parameter.  The  MLM  routine  was  run  five  times  on  a  two  second  set  of  data  that 
included  the  IP  arrival  at  8  Hz  in  Figure  5-4a.  For  each  run,  ^  was  changed, 

A 

its  value  ranging  from  .01  to  .05  in  increments  of  .01.  The  results,  10  log  P, 

plotted  in  figure  6-7,  varied  by  5  db,  due  to  bias  dependence  on  V  .  With 

* 

bias  corrections,  the  corrected  levels,  10  log  Pg^^  ,  remain  within  1  db  of 
each  other. 


Source  Level 

For  each  shot  point  in  the  Rose  experiment,  the  explosive  weight  and 
an  estimated  shot  depth,  based  on  sinking  rate  and  source  monitor  times,  were 
known.  Employing  this  information  in  an  empirical  relation  (Wakely,  1977),  an 
expression  for  the  pressure  waveform  at  a  range  R  was  computed.  The  model 
waveform  includes  four  bubble  pulse  periods  following  an  initial  shock  pulse. 
The  waveform  was  Fourier  transformed  and  the  squared  magnitude  of  the 
resulting  components  were  averaged  across  the  estimation  bands  in  the  data. 


25  pound  charge  detonated  at  35,9  meters  obtained  with  Wakely  formula 
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After  correcting  back  to  a  1  meter  range,  the  results  became  the  source 
pressure  levels  used  in  the  TL  calculations.  In  Fig.  6-5,  representative 
temporal  and  spectral  profiles  for  a  25  lb  charge  detonated  at  35.9  meters  are 
shown.  For  this  shot,  the  averaging  in  frequency  was  particularly  important  in 
the  8  Hz  band  because  of  a  10  db  dip  in  the  level  near  7.5  Hz.  Since  the  shock 
wave  intensity  in  all  shots  used  in  our  work  was  above  the  cavitation  limit, 
the  energy  near  the  sources  was  incoherent  and  a  Lloyd  mirror  correction  at 
the  shot  points  would  not  be  valid.  Source  level  behavior  from  a  cavitating 
shock  wave  at  the  surface  is  a  nonlinear  acoustics  problem  which  still  needs 
analysis. 

Implemented  transmission  loss  equation 

Applying  the  five  corrections  outlined  above,  the  expression  for 
partitioned  transmission  loss  In  eq.  6-11  can  be  rewritten  as: 

TL^4'l=S/-(p-[(oAget('(<'t)-2Oi^S./C)i^e-2.o^L^"2oAgA!:0 

(6-16) 

where  S  denotes  hydrophone  sensitivity,  e  is  the  sampling  correction,  and  LM 
and  AD  are  the  mirror  effect  and  array  directivity  corrections,  respectively. 
LM  and  AD  effects  were  not  applicable  to  the  vertical  array  (MABS)  estimates. 

A 

P  (f,  v')  represents  the  magnitude  of  the  frequency -wavenumber  estimate  after 
bias  correction.  Results  obtained  from  the  use  of  this  expression  are  now 
presented. 
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iransmission  Loss  Lstimates  for  IP  events 

In  figure  6-6,  TL  estimates  at  8  Hz  for  the  first  major  event  (IP)  in 

both  MABS  and  ESP  data  from  lines  2S  and  2L  are  shown  up  to  an  offset  of  42 

kilometers.  Line  2S  shots  consisted  of  alternating  5  and  25  lb  charges.  Only 
25  lb  shots  are  used  in  the  processing  of  MABS  data  from  this  line,  but  both 
sizes  are  included  in  the  estimates  for  the  ESP  data,  with  the  5  lb  shots 

being  indicated  in  the  figure.  Line  2L  at  these  offsets  used  180  lb  charges 

exclusively.  The  following  points  should  be  noted: 

i)  Line  2S  shots  with  5  lb  charges  generally  have  higher  TL  estimates 
than  those  using  25  lbs.  A  bolder  line  is  drawn  through  points  corresponding 
to  2b  lb  charges  only.  The  approximate  5  db  difference  in  the  calculated  TL 
for  5  and  25  lb  shots  is  probably  due  to  errors  in  source  level  estimates.  The 
depths  at  which  both  size  charges  detonated  were  roughly  the  same  (40  m.),  so 
tnat  Lloyd  mirror  corrections  on  the  source  levels,  even  if  they  were  valid, 
would  not  change  the  relative  difference  in  the  estimates. 

ii)  If  we  ignore  the  5  lb  data  and  note  that  line  2L  was  processed 
with  24  channels,  the  ESP  2S  data  with  12,  and  the  MABS  with  8,  we  can  see  a 
steady  progression  in  the  TL  curves  with  the  data  with  the  least  number  of 
sensors  having  the  least  loss.  The  difference  in  the  estimates  for  line  2L 
data  and  2b  data  could  again  be  due  to  error  in  source  level  estimates,  line 
2L  having  used  180  lb  charges.  The  approximate  5db  mean  difference  of  the  25 
Id.  MrtbS  and  ESP  2S  data,  however,  cannot  be  caused  by  the  use  of  erroneous 
source  levels.  The  discrepancy  may  be  attributed  to:  a)  the  fact  that, 
geographically,  the  MABS  and  ESP  2S  data  are  samples  of  different  locations 
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Results  of  transmission  loss  calculations  on  first  arrivals  (IP) 
at  8  Hz  for  ESP  lines  21  and  2S,  Including  vertical  (MABS)  array  data 
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and  the  crust  is  not  laterally  homogeneous;  b)  statistical  fluctuations  in  the 
MLM  estimator;  and  c)  errors  in  MLM  bias  correction.  Since  the  bias 
corrections  employed  are  based  on  empirical  results  and  since  MLM  bias  is  very 
sensitive  to  the  number  of  sensors  used,  the  latter  is  probably  more 
significant.  If  we  look  at  fig.  6-7,  however,  in  which  TL  estimates  for  the 
same  set  of  data  are  shown  without  bias  corrections,  we  can  see  that  the 
variation  of  the  estimates  with  the  number  of  sensors  used  has  been 
significantly  reduced  in  the  corrected  set. 

iii)  There  is  a  consistent  drop  in  TL  for  all  3  data  sets  between  the 
ranges  of  25  and  40  km.  As  we  have  mentioned,  this  is  often  encountered  in 
refraction  data.  This  drop  is  about  6  db  In  2S  ESP  data,  10  db  in  2S  MABS,  but 
only  2  db  in  the  2L  data.  The  latter,  however.  Is  undersampled  so  that 
evidence  for  a  greater  focusing  effect  between  offsets  of  33  to  40  km  may  have 
been  missed. 

Figure  6-9  illustrates  the  results  produced  when  TL  estimates  for  line 
2L,  out  to  a  range  of  104  km,  are  calculated  for  4  separate  frequency  bands 
and  are  "corrected"  for  geometrical  spreading.  A  value,  10  log  rl(r  being  the 
horizontal  offset  in  meters),  was  subtracted  from  each  TL  value.  In  this 
drawing,  an  ideal  pressure  wavefront  with  a  simple  spherical  attenuation  would 
appear  as  a  horizontal  line.  We  observe  that  the  actual  loss  in  the  crust 
increased  with  range  somewhat  faster  than  the  rx  dependence.  Assuming  the 
geometrical  factor  has  been  accounted  for,  this  added  loss  reflects  the 
absorption  of  energy  that  has  taken  place  along  the  path.  In  this  figure,  the 
"resonance"  phenomenon  of  a  low  point  in  TL  between  offsets  of  30  and  40  km 


TRANSMISSION  LOSS  (d  b ) 
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appears  more  dramatically  for  line  2L  data  than  in  figure  6-7  because  of  the 
geometrical  spreading  factors  applied. 


Parameters  of  linear  regression  applied  to  TL  estimates 

In  Table  6-2,  the  results  of  applying  a  linear  regression  to 
transmission  loss  estimates  versus  offset  are  presented.  The  parameters  of  the 
regression  tabulated  are: 


N 


the  number  of  TL  estimates  used  In  the  regression; 


(db/km)  the  slope  of  the  fitted  line  with  respect  to  offset, 

indicating  loss  above  (or  below)  geometrical  losses  due 
to  absorption  and  other  effects; 

G>  (db)  the  intercept  at  zero  offset  of  the  regression  which  is 
an  indication  of  "Insertion  losses"  such  as  reflection 
losses  at  layer  boundaries  and  transmission 
coefficients; 

X  the  "coefficient  of  determination"  or  correlation 

coefficient  indicating  the  quality  of  fit  achieved  by 
the  regression.  Values  closer  to  1  Indicate  a  better 
fit  than  values  near  zero; 


6 


the  standard  deviation  of  the  regression 


Q 


a  dimensionless  attenuation  factor:  the  ratio  of  energy 
stored  in  one  cycle  to  the  energy  lost  during  the 

20i^e  (Clay  *  Medwin,  1977) 


cycle: 


•  I  •  W..V  VI  V  I 


Regression  results  for  path  IP  from  line  2L  in  table  6-2  are  presented 
for  the  entire  line  and  also  separately:  i)  for  offsets  up  to  35  km;  and  1i) 
for  offsets  beyond  35  km.  The  correlation  coefficients  for  the  latter  two  sets 
are  consistently  higher  than  for  estimates  made  from  the  entire  line,  since 


Data  Set 

Freq.  (Hz) 

H 

(db/km) 

ft  (db) 

-e- 

Line  2L 

5 

11 

.12 

34.59 

.53 

3.22 

(20  to 

8 

11 

.17 

38.13 

.66 

3.38 

104  km 

11 

10 

.09 

45.29 

.54 

2.01 

offset) 

14 

10 

.11 

50.66 

.63 

2.06 

Line  2L 

5 

3 

.45 

28.33 

.98 

.34 

(20  to 

8 

3 

.22 

38.9 

.94 

.3 

35  km 

11 

3 

-.08 

50.8 

1.00 

.13 

offset) 

14 

3 

-.13 

57.5 

.98 

.16 

Line  2L 

5 

8 

.20 

28.56 

.69 

2.85 

(35  to 

8 

8 

.24 

32.56 

.7 

3.39 

104  km 

11 

7 

.15 

40.58 

.73 

1.77 

offset) 

14 

7 

.15 

47.54 

.65 

2.12 

Line  2S 

8 

6 

-.13 

44.43 

.08 

2.11 

(ESP  5# 

12 

6 

-.93 

61.05 

.99 

.53 

shots;  16 
to  30  km) 

16 

6 

-.67 

60.48 

.65 

2.35 

Line  2S 

8 

3 

.23 

34.90 

.72 

.81 

(ESP  25# 

12 

3 

.18 

37.08 

.08 

3.39 

shots;  17 
to  30  km) 

16 

3 

-.50 

58.50 

.98 

.45 

Line  2S 

8 

11 

.77 

22.06 

.78 

2.23 

(MABS;  8 

12 

11 

-.02 

40.36 

.03 

2.63 

to  25  km 
offset) 

16 

11 

-.23 

47.46 

.37 

1.58 

TABLE  6-2 


Results  of  linear  regression  for  path  IP 
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the  effect  of  the  nonlinear  behavior  at  "resonance"  near  35  km  is  decreased. 
This  division  of  2L  data  Is  also  motivated  by  the  fact  that  data  below  35  km 
may  reflect  losses  in  propagation  through  "layer  3",  while  estimates  made  from 
data  beyond  35  km  involve  energy  that  has  Interacted  with  the  Moho. 

The  most  consistent  feature  in  the  results  for  2L  in  table  6-2  is  the 
increase  in  insertion  loss,  (3  ,  with  center  frequency.  This  loss  Increases 
an  average  of  9  db  for  each  3  Hz  Increment  in  frequency.  The  magnitude  and 
frequency  behavior  of  @  is  essentially  on  agreement  with  the  results  of 
Baggeroer  and  Falconer  (1981)  for  estimates  of  transmission  loss  for  events 
interpreted  as  Moho  refractions  in  the  CANBARX  experiment.  The  slopes,  c<  , 
however,  averaging  about  .2  db/km  for  line  2L  "Moho"  data,  are  consistently 
lower  than  those  in  the  CANBARX  paper  (which  is  on  the  order  of  .5  db/km)  and 
are  also  lower  than  attenuations  extrapolated  from  data  published  by  Hamilton 
(1972),  which  tends  to  be  closer  to  1  db/km.  The  data  in  the  latter  paper  is 
relatively  sparse  at  low  frequencies  and  concerns  propagation  in  marine 
sediments.  Acoustic  behavior  in  basement  basalt  and/or  Moho  would  be  expected 
to  be  more  efficient  as  the  present  results  Indicate. 

As  with  the  Rose  data,  the  estimates  in  the  CANBARX  paper  are  based  on 
crustal  refraction  data  and  the  procedure  used  in  evaluationg  TL  is  the  same 
as  used  here  except  that  MLM  bias  corrections  were  not  performed.  The  fact 
that  large  bias  corrections,  essentially  based  on  empirical  results,  were  used 
in  the  24  channel  2L  data  and  that  similar  corrections  (although  they  would  be 
smaller  for  the  CANBARX  array,  with  a  maximum  of  12  functioning  channels)  were 
not  applied  to  the  CANBARX  data  may  account  for  the  discrepancy  between  the 
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two  estimates  of  c<,  .  CANBARX  results  are  also  based  on  data  taken  from  five 
experiments  while  the  ROSE  2L  estimates  involve  data  from  up  to  eleven 
separate  shots. 

The  values  of  Q  obtained  here,  in  keeping  with  lower  estimates  of 
in  ROSE  results,  are  higher  than  In  CANBARX.  As  such,  they  tend  to  be  closer 
to  the  results  of  Jacobson  (Jacobson  et  al,  1981)  in  which  values  of  Q  '  ,  the 
“specific  quality  factor",  obtained  in  a  sedimentary  region  in  the  Bay  of 
Bengal,  approached  values  less  than  .01  at  the  greatest  depths  sampled. 
Estimates  of  Q  are  lacking  in  the  tables  for  data  in  which  the  resulting  c<. 
estimate  was  negative.  These  negative  estimates  occur  In  lines  in  which  the 
attenuation  of  energy  was  less  severe  than  that  caused  by  spherical  spreading, 
possibly  due  to  the  effect  of  "resonance". 

ROSE  line  2L  data,  although  relatively  sparse  with  respect  to  shot 
density,  is  considered  to  be  of  higher  quality  than  the  line  2S  data  in  which 
a  smaller  number  of  data  channels  was  used.  In  table  6-2,  results  for  path  IP 
are  also  presented  for  ESP  and  MABS  data  for  2S.  The  parameters  of  the 
regression  are  considerably  more  scattered  than  those  for  2L  and  the 
correlation  coefficients  of  some  2S  results  decreases  below  .1.  In  comparing 
the  2L  and  2S  parameters,  with  the  data  from  2S  involving  offsets  up  to,  but 
not  including,  the  ranges  at  which  "resonance"  occurred,  a  consistent  feature 
appears  to  be  the  decrease  of  the  estimated  slope,  ,  with  increasing 
center  frequency.  In  both  the  25  lb  ESP  and  the  MABS  2S  cases,  o<  is  a 
maximum  at  the  lowest  frequency  and  decreases  so  that  at  the  upper  center 
frequency,  it  is  negative.  This  Is  a  suspect  result  since  ,  reflecting 
(at  least  partially)  absorption  losses  in  the  crust,  would  be  expected  to 
Increase  with  frequency,  perhaps  linearly  as  in  the  data  published  in  Hamilton 
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Path 

Data  Set 

Freq.  (Hz) 

N 

o<  (db/km) 

g(db) 

$ 

_Q_ 

2P 

Line  2L 

5 

11 

.04 

r 

37.91 

\ 

.11 

2.92 

487 

(20  to 

8 

11 

.07 

41.25 

.26 

3.11 

445 

104  km 

11 

10 

.04 

46.53 

.3 

1.57 

1072 

offset) 

14 

10 

.08 

49.62 

.34 

2.73 

682 

2P 

Line  2S 

8 

5 

-.20 

45.74 

.06 

3.06 

(ESP  5# 

12 

5 

-.45 

50.58 

.46 

1.92 

- 

shots;  16 
to  30  km) 

16 

5 

-.26 

51.06 

.29 

1.59 

* 

2P 

Line  2S 

8 

3 

.67 

21.76 

.32 

3.43 

46 

(ESP  25# 

12 

3 

.83 

16.84 

.99 

.21 

56 

shots;  17 
to  35  km) 

16 

3 

.28 

39.36 

.18 

2.06 

222 

2P 

Line  2S 

8 

8 

.59 

24.18 

.4 

2.43 

53 

(MABS;  16 

12 

8 

-.29 

44.00 

.22 

1.66 

- 

to  25  km) 
offset) 

16 

8 

.19 

35.56 

.03 

2.99 

328 

1W 

Line  2L 

8 

5 

-.005 

28.75 

.0008 

1.93 

(20  to 

11 

5 

.26 

32.24 

.26 

2.16 

769 

52  km) 

14 

5 

.87 

47.59 

.87 

2.12 

292 

1W 

Line  2S 

8 

12 

.1 

29.9 

.39 

.95 

1455 

(ESP;  15 

12 

12 

-.07 

29.75 

.04 

2.19 

- 

to  42  km; 
all  shots) 

16 

12 

.07 

19.34 

.02 

.98 

4158 

1W 

Line  2S 

8 

6 

.05 

31.87 

.19 

.8 

2910 

(ESP  to  42 

12 

6 

-.19 

34.6 

.28 

2.13 

- 

km;  25#  shots) 

1U 

Line  2S 

8 

6 

.10 

29.48 

.22 

.89 

1455 

(to  30  km; 
5#  shots) 

12 

6 

-.16 

30.74 

.3 

1.20 

- 

"Layer 

Line  2S 

8 

4 

-.1 

54.34 

.13 

.93 

2" 

(ESP) 

12 

4 

-.37 

58.71 

.44 

1.59 

- 

16 

2 

-.58 

67.5 

1.00 

- 

Table  6-3 

Results  of  linear  regression  for  paths  2P,  1W,  and  "layer  2"  returns 
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(1972).  Although  the  2S  data  is  scattered,  this  pattern  is  corroborated  in  the 
2L  results  within  the  20  to  35  km  offset  range.  This  phenomenon  may  be  caused 
by  the  "resonance"  affecting  the  estimated  slope  at  offsets  less  than  35km. 
Another  possibility  is  that,  as  a  general  rule,  bias  corrections  tend  to  be 
largest  at  higher  frequencies  where  original  frequency/wavenumber  estimates 
are  usually  lower.  Higher  frequency  data  may  have  tended  to  be 
“over-corrected",  causing  this  pattern  in  the  regression  parameters.  The 
pattern,  however,  seems  to  be  associated  with  data  collected  up  to  the 
resonance  point  only.  Line  2L  data,  beyond  35  km,  is  not  affected  as  severely. 

In  table  6-3,  results  of  regressions  performed  on  TL  estimates  for 
paths  2P,  1W,  and  the  low  level,  possibly  "layer  2",  events  are  tabulated.  In 
accord  with  results  in  the  CANBARX  paper,  the  insertion  losses  for  the 
multiple  reflection/refraction  path,  2 P,  are  greater  than  those  for  the 
primary  IP,  but  attenuations,  ,  especially  In  the  24  channel  data,  tend 
to  be  smaller  (near  zero).  The  phenomenon  of  relatively  high  multiple 
amplitudes  is  commonly  observed  in  refraction  profiles  but  remains  to  be 
analysed  rigorously.  In  all  the  2P  data  presented,  the  correlation 
coefficients  are  considerably  smaller  than  those  for  IP,  possibly  due  to  the 
rough  topography  encountered  in  the  waterborne  segments  of  these  events. 

The  effects  of  seafloor  topography  also  extend  to  the  1W  data  in  which 
all  of  the  regression  parameters  tabulated  again  have  relatively  low 
correlation  coefficients.  As  discussed  previously,  the  theoretical  path  for 
these  events  encounters  the  seafloor  before  arriving  at  the  receivers.  For 
line  2L  data,  insertion  loss  and,  in  this  case,  increases  with  frequency. 
Care  must  be  taken  in  interpreting  this  result  because,  as  we  have  seen,  large 
Lloyd  mirror  and  directivity  corrections  are  incorporated  in  the  higher 
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frequency  estimates.  Although  the  quality  of  the  linear  fit  is  poor,  due  to 
the  rough  bathymetry,  the  estimated  slopes  from  line  2S  average  near  zero, 
implying  near  spherical  spreading  for  1W. 

Finally,  although  only  four  events  could  be  identified  as  being 
possible  "layer  2"  arrivals,  the  regression  results  for  theses  available  only 
for  line  2S,  show  the  same  decreasing  slope  pattern  with  frequency  that  was 
encountered  in  Table  6-2.  The  magnitudes  of  the  slopes  roughly  correspond  to 
the  slopes  obtained  for  IP  in  ESP-2S,  but  the  insertion  loss  is  markedly 
higher.  This  is  in  sharp  contrast  with  the  low  <**.  for  events  considered  to  be 
layer  2  arrivals  in  the  CANBARX  paper.  The  identification  of  this  set  of 
events  in  the  ROSE  data  remains  undetermined:  the  arrival  times  correspond!’ ng 
to  those  of  possible  shear  arrivals  (which  would  account  for  the  high 
insertion  loss  due  to  poor  coupling  between  congressional  and  shear  waves  at 
the  basement),  but  the  estimated  phase  velocities  at  the  array  tending  to  be 
too  high  for  shear  propagation. 

Summary 

In  applying  MLM  array  processing  techniques  to  the  analysis  of  data 
from  the  ROSE  experiment,  the  effectiveness  of  the  velocity  estimation 
procedure  was  diminished  by  extremely  rough  topography.  Estimated  phase 
velocities  at  the  array  reflect  variations  In  bathymetry  as  well  as  the 
properties  of  crustal  sections  with  which  generated  energy  interacts.  Results 
obtained  in  applying  the  same  analysis  techniques  to  data  from  experiments 
such  as  CANBARX  (Baggeroer  &  Falconer,  1981)  and  FRAM  II  (Duckw  th  et  al, 


Iy82),  in  which  bottom  roughness  was  much  less  severe,  have  been  more 
successful  in  the  inversion  of  the  velocity  information  obtained  into  higher 
resolution  crustal  models.  We  have  shown  that,  even  in  a  "worst  case" 
situation,  the  ability  of  the  algorithm  to  discriminate  arrivals  by  means  of 
relative  phase  velocity  estimates  is  valuable  when  combined  with  more 
conventional  methods  of  processing  refraction  data. 

Apart  from  considerations  of  velocity  structure,  we  have  shown  that 
amplitude  information  obtained  from  MLM  estimates  can  be  used  effectively  for 
obtaining  estimates  of  transmission  loss  in  the  crust.  Although  some  work  has 
been  done  in  determining  TL  in  marine  sediments,  the  work  described  here  and 
in  the  CANBARX  paper  (Baggeroer  &  Falconer)  is  a  rare  attempt  at  estimating 
crustal  losses.  The  results  of  both  papers  agree  In  Insertion  loss  estimates 
for  primary  crustal  paths.  The  attenuation  of  primary  paths,  on  the  order  of 
.5  db/km  in  the  CANBARX  results,  is  somewhat  smaller,  on  the  order  of  .2 
db/km,  in  the  work  done  with  the  ROSE  data. 

The  ROSE  results  incorporate  corrections  obtained  by  the  introduction 
of  a  procedure  for  the  removal  of  bias  effects  In  the  MLM  algorithm,  so  that 
TL  estimates  are  more  accurate.  This  procedure,  based  on  empirical  results,  is 
applicable  to  MLM  estimates  obtained  from  sparse  arrays,  which  are  often  the 
only  economical  and  practical  means  of  obtaining  the  benefits  of  array 
techniques  in  the  marine  environment. 
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The  correlation  function  of  a  simple  random  process  x(t)  is  defined  as: 

R*  Ct)  -  £  £ octe)  U -t)'] 

where  denotes  expectation,  and  the  superscript  *  is  the  complex  conjugate, 
in  the  situation  shown  in  fig.  4-1,  a  process  x(t)  is  the  input  to  a  pair  of 
linear  filters  with  frequency  responses  H,(f)  and  Ht(f),  and  impulse  responses 
h.lt)  and  ha(t).  The  output  of  a  filter  is  just: 

ajl(£)  ~  oc  C  £)  *  JsJ.it) 

-  J"  ocCV)  JJ(Jz-'r)  At- 

where  *  denotes  convolution.  The  correlation  function  of  the  output  process 
yt(t)  with  any  desired  process,  d(t),  is: 

Os)  *  £  { 

—  «D 

c  j  JLX  (k"r)  Ct-<&+6)  c&n" 

With  the  substitution  'Y>  =  /fc-‘T'  ,  we  get: 

(e)  -  j  J^,  6r')  (e-'r^Ar1 

v  .dO 

~  (A4-1-1) 

Likewise,  the  correlation  function  (  6  )  is: 


R G)  -  *(*-«)} 

-  j  X* (t-t- r)  € £<SC* W* (t)}  <4t 

With  the  substitution  /y-1  =  »  we  have: 

R&i'Ct)  -  j - s)  R<^Ct')At 

*  JL,  (-s')  * 


(A4-1-2) 
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If  we  take  the  Fourier  transform  of  A4-1-1  and  A4-1-2,  we  obtain  the 
spectral  density  functions: 


From  A4-1-3,  using  d(t)  =  y,(t),  we  obtain: 


(A4-1-3) 

(A4-1-4) 


S*. ({)  -KCQS^U) 

jsing  d(t)  =  x(t)  in  A4-1-4,  we  get: 

(V  -  aqM'W&q)- iwi'&il.,. 

,  if  d(t)  is  y^(t)  in  A4-1-3: 

Swxty*  /// (p ({) 

letting  x(t)  be  d(t)  in  A4-1-4: 


(A4-1-6) 
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wnere: 


Appendix  4-2 

calculation  of  spectral  density  function  of  the  output  of  an  array  processor 
Referring  to  fig.  4.2,  the  output  of  the  processor  is  just: 

y.Ct)  - 

The  correlation  function  of  this  output  is  then: 

Ry.  61  *  £ 

Taking  the  expectation  inside  of  the  summations,  we  obtain: 

=  ?  f 

we  rewrite  equations  A4-1-1  and  A4-1-2  from  Appendix  4-1: 

(s')  =  3-»  ^  to 

Rju,,  6>)  *  3^-6)*$^  fr) 

^  n*t- c-o/ 

where  dU)  is  any  desired  process.  Identifying  y. (t)  with  d(t)  in  A4-2-2, 
equation  A4-2-1  becomes:  £ 

a) 

how  identifying  d(t)  with  x. (t),  we  obtain: 

Rf  <a)  =  f  ^  5. to  *  5,  to  *  £***,•  6) 

t.i.  4-u^  r 4. . 
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Taking  the  Fourier  transform,  we  get: 
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Uerivation  of  Optimal  Response  Function  for  Arbitrary  Noise  Covariance 
Vie  wish  to  minimise  with  the 

constraint  that  &({.**)  •  '  -  %  • 

Via  Lagrangian  Multipliers: 


MIN:  S 


b.lf) 


&J{)Si  qK*(p,x({){i,  <^)e *  - 1} 


a)  Take  partials  >  and  set  to  zero: 

. >.  «*m««m** 


-  a  I'm*  A  Ac. 


=  0 


b)  Define  vector  6j[f)  =  L 


-  ie 


jlTr* a, 


]»  the  steering  vector 


£  ],  and  the  KxK  Hermitian  matrix 


LS-  (f)J,  so  that  we  may  write  the  above  as 

_  j  a  \  1  y"  X  /  A  \  ^ " 


*x tSj  ^r|  <2.  *■ \q)  **) "  ® 

or  -I 

(zty)  -  -  ISp  ({U  £*({•**)  (by  taki"9  conju^transpose) 

c)  Take  partial  -hr®  and  set  to  zero: 

k &xq) „  qt£ *  =± 

d)  Substitute  for  £  from  A4-3-10: 

±  -  - 

e)  Now  substitute  -  — £  into  equations  A4-3-1: 


(A4-3-2a) 


f)  bince 
we  have: 


<k  =  ^({a)  ' £*({**> 
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q]  The  spectral  density  function  of  the  output,  which  is  also  the  MLM  estimate 
of  the  frequency  wavenumber  function,  is  found  by  substituting  from  eqs. 

A4-3-2  into  the  quadratic  form: 


^q)-s.-rtsiqa^‘ 
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Appendix  4-4 

uerivation  of  Optimal  Response  Function  for  Uncorrelated  Sensor  Noise 


we  wish  to  minimize  S^(f)  =  ^ I  SttfllAtWl2- 
constraint  that  =  1  =  Jl  Gity)  e 

Via  Lagrangian  multipliers: 

£*({)£'  l] 


with  the 


MIN: 

over 

^(f) 


(A4-4-1) 


(A4-4-2) 

(A4-4-3) 


—  - - — 

a)  Take  partials  and  set  to  zero: 

a  SZq)  <£0+V#  *fc-0 

or  r  -Mpe.**2™***- 

^JkSv  -  sX,  ® 

b)  Take  partial  /dAi^)  and  set  to  zero: 

£  6.^)  1  -o 

c)  Substitute  value  of  bj^f)  in  A4-4-2  into  A4-4-3,  and  obtain 

2?r  t  z<^>  /  x  u 
Xf^)  =  -  a 

d)  Substitute  value  of  X(f)  into  A4-4-2  and  get  final  results: 
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Appendix  4-5 

tvaiuation  of  correlation  function  of  Fourier  coefficients  from  a  sensor  pair 

Let  X  If  )  =  (  oc .x(£)aa>C*) 

x~>  1 


be  the  Fourier  coefficient  at  f,  ,  associated  with  the  windowed  time  series 
from  sensor  i.  The  window  lenath  is  T  sec.  Then  thp  correlation  function* 


Substitute  the  density  function  associated  with  (t,  ,tx  ): 

» So  **<**)  **■ 

» J A*r W 

where  l*(f)  = 


if  tiw  is  the  effective  bandwidth  of  the  transform  of  the  window  function 
(which,  for  a  cosine  window  is  8/3T  'fts  2.6  Hz  for  T  =  1  sec),  then: 


r 


o 


(f,  -fz  >  BW) 
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Appendix  6-1 

Corrections  to  Frequency/Wavenumber  Estimates  due  to  Sampling 


For  a  pressure  waveform  at  location  r  of  length  T,  with  an 
approximately  constant  spectral  density,)*^?,  r), 
can  write:  * 


across  frequency  band,  W,  we 


/V  i. 


=  T  (*±}(jt))  *  VjJf  (Jo,  4) 

—  T/t. 

or: 

a  TfL 

7*/  ~  Wt 

if  we  sample  data  at  an  interval  aT,  we  can  approximate  this  as: 

A  /  jll*/  $J 

(4,^aT)«.(»  Cut  ^,-^.up 

where  N  =  T/aT.  In  terms  of  the  two  sided  discrete  Fourier  transform  of  the 
raw,  broadband  mean  square  pressure,  n  T),  this  becomes: 


£L  Of1 

UJT  $T0 


•  2-rr^iit  *  z' 


where  k ,,  is  the  frequency  term  number  corresponding  to  f„  ,  and  af  is  the 
coefficient  spacing  in  Hz.  Rewriting  this  as: 

and  recognizing  that  the  quantity  in  brackets  is  zero  unles  k  =  iN 
ti  integer),  we  have: 
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For  data  sampled  exactly  at  the  Nyquist  rate,  1/aT  =  2W  ,  where  W  is  the 
total  bandwidth  of  the  input,  the  spacing  of  the  N  (=  T/aT)  frequency 
coefficients  in  the  OFT  is  just: 

_  avu  Vr _ !_ 

NJ  T  "■  T 

The  number  of  terms  (M)  in  the  present  band,  W,  is  then: 

M  =  "//T  =  WJT 

we  can  rewrite  the  expressions  for  the  density  function  as: 
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jendix  6-2 


Lloyd  Mirror  Correction 


deterring  to  tig.  fa-1,  EDA  represents  a  pure  sinusoidal  plane  wavefront  with 
frequency  f,  arriving  at  a  sensor  at  point  A  at  time  t.  The  vertical  angle 
is  oc  and  hydrophone  depth  is  BA  =  d.  The  sensor  will  also  be  influenced  at 
tnis  moment  by  a  surface  reflected  arrival  that  has  traversed  the  extra 
distance  L)C  +  CA.  Assuming  specular  reflection: 


=  (3, 


tne  geometry  of  the  situation  will  be  as  shown.  The  extra  distance  travelled 
is  then: 

DC  +CA  - 


~  c<. 

it  the  upcoming  arrival  is  a  pure  sinusoidal  plane  wave  with  frequency  f  and 
unit  amplitude,  the  waveform  at  A  can  be  written  as: 

ico^r  - 


n 


where  we  have  assumed  perfect  reflection  at  the  surface  except  for  a  180‘ 
phase  shift. The  amplitude  of  the  resultant  waveform  is  then: 


/-  e 


1  dAstoC 
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12  Defense  Technical  Information  Center 
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ATTN :  DCA 
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Naval  Oceanographic  Office 
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Bay  St.  Louis,  MS  39522 
1  ATTN:  Code  8100 

1  ATTN:  Code  6000 

1  ATTN:  Code  3300 

1  NODC/NOAA 
Code  D781 

Wisconsin  Avenue,  N.W. 

Washington,  D.C.  20235 

1  Mr.  Michael  H.  Kelly 

Administrative  Contracting  Officer 
Department  of  the  Navy 
Office  of  Naval  Research 
Eastem/Central  Regional  Office 
Building  114,  Section  D 
666  Summer  Street 
Boston,  MA  02210 
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